texts: Black & White / Greyscale Maps and Illustrations

texts:
I think greyscale maps & illustrations (with geographic data) seem to be abandoned and the colour maps are mostly used nowadays. This is understandable and personally I also need less time to produce a colour map than a grey scale map, although in the end I fell much more satisfied with the grey scale results.

What are your thoughts?
Is it still feasable to elaborate grey scale maps and illustrations?
Where can I find some good tutorials on this(e.g. advises on grey color selection, how to deal with such a limited color range,etc)?

I am not looking for "converters" from color to grey scale, like the one from ESRI (that seems to work ok) that can easily and with soem quality export a colour map to grey scale: http://mappingcenter.esri.com/index.cfm?fa=ask.answers&q=2163


texts: Color different in Maps exported from ArcMap to PDF and other formats?

texts:
When exporting the same map in ArcMap into different formats the color is not consistent.

The map is a roadway network and I use grey 30% solid line as the symbol in ArcMap.

However, when i export it to PDF and all the other formats, only the PDF resembles the color in ArcMap. All the rest appear like a color of grey 10%.

My problem is I have to generate PDF and picture format (which can be directly inserted into Word document) maps at the same time, but the colors appear different.

Instead of spending time to make two different map document with adjusted colors, what can I do?

The first image below is a jpeg format map, you can hardly see the grey color roadways in the background. The second one is the pdf format map, the grey color shows up exactly as the same as what I see in the ArcGIS.


texts: Classifying small structure in water using IR band

texts:
I was wondering how to identify small structures (eg. beacons) in a water body? I have worked with RGB images and its almost impossible to classify waterbodies in these images.

How does the IR band help here? Will they be good enough to say that this pixel belongs to water body and this pixel doesn't?

I hope I am clear enough. Kindly ask if you need more information.


texts:
This is my first time writing in this forum and I hope this question offers some quality of information easing a solution. Here goes.

I am new in analyzing hydrology with DEMs and I have followed several steps to identify watersheds and drainage basins. (I am using QGIS for all my operations)

Obtained the DEMs and combined them using "Build Virtual Raster"
Filled the holes using "Fill sinks (wang & liu)" providing a filled
DEM and Flow Directions
With the filled DEM I've calculated the Strahler order (which I'm not familiar with) and identified rivers with this tool.
A new raster was calculated with the Strahler order's values >= 8 resulting in an acceptable river raster layer (validated through a basemap).
Then I've done a Channel network and drainage basins with a threshold 8 (to have the same acceptable result as my Strahler order) and saved the channels and drainage basins. (I suspect my error might be in this step)
Then the result that I get gives me a drainage basin going straight through a lake, which makes no sense I find. Shouldn't the lake it'self be a receiver of a whole network of watershed?

Here is a screenshot of the result I'm getting with the DEM in the background.

The drainage basins are the red limits. I figured the limits of the drainage basins would either be where I inserted the arrows 1 or 2 or both, as we can clearly see the lake/reservoir in the DEM. I figure I might have missed a step to tell my procedure the limit of every water plan in this process as it doesn't seem to take them into account.


texts:
I want to filter lines which represent ship movement in a manner, that linestrings where the distance between two vertexes is greater than 1000 km are removed.
I have a piece of code, however I am getting the error message that the function ST_Distance does not exist.

> ERROR:  function st_distance(geometry_dump, geometry_dump) does not exist
LINE 1: SELECT ST_Distance(ST_DumpPoints(a.line), ST_DumpPoints(b.li...
               ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
SQL state: 42883
Character: 8

Here is my code:

SELECT ST_Distance(ST_DumpPoints(a.line), ST_DumpPoints(b.line)) dist
FROM ais_2017_09_move a, ais_2017_09_move b
WHERE mmsi = 566035000
AND dist < 1000;


texts:
I would like to produce series of maps with three different base maps. Within each data frame I have the same point shapefile with 1000 locations. Data driven pages allows only to produce maps from one frame. This means that in the resulting pdf only locations from first frame is changing and not for other two frames. I found that this could be done using ArcPy.

Could someone help me please with a simple example?

Are there also another possibilities without coding?


texts: ArcPy UpdateCursor deleting duplicate records?

texts:
I am working with .dbf files. I have two columns of data that have duplicates. A column of parcel numbers with duplicates and a column of letters that are duplicated. I want to delete all parcel numbers first and then delete all duplicate letter numbers after. Right now, I'm just repeating code twice.

Is there an easier way?

I'm brand new at this so that's why the code is repeated and looks so chunky. It's code that I mostly found off this site previously.

Here's what I'm working with:

import arcpy
import os
from arcpy import env

inShapefile = "M:/GIS/LEA Maps/2019/List/GISLEA.dbf"
checkField = 'PARCELNO'
updateField = 'dplicate'

def findDupes(inShapefile, checkField, updateField):
    with arcpy.da.SearchCursor(inShapefile, [checkField]) as rows:
        values = [r[0] for r in rows]

    with arcpy.da.UpdateCursor(inShapefile, [checkField, updateField]) as rows:
        for row in rows:
            if values.count(row[0]) > 1:
                row[1] = 'Y'
            else:
                row[1] = 'N'
            rows.updateRow(row)

if __name__ == '__main__':
    fc = r'M:/GM/LEA/LEADR.dbf'

    findDupes(fc, checkField, updateField)

arcpy.DeleteIdentical_management(fc, 'PARCELNO')

inShapefile = 'M:/GM/LEA/LEADR.dbf'
checkField = 'LetterOrde'
updateField = 'dplicate'

def findDupes(inShapefile, checkField, updateField):
    with arcpy.da.SearchCursor(inShapefile, [checkField]) as rows:
        values = [r[0] for r in rows]

    with arcpy.da.UpdateCursor(inShapefile, [checkField, updateField]) as rows:
        for row in rows:
            if values.count(row[0]) > 1:
                row[1] = 'Y'
            else:
                row[1] = 'N'
            rows.updateRow(row)

if __name__ == '__main__':
    fc = 'M:/GM/LEA/LEADRlo.dbf'

    findDupes(fc, checkField, updateField)

arcpy.DeleteIdentical_management(fc, 'LetterOrde')


texts: How to add additional HTTP response headers?

texts:
I installed recent MapProxy (1.7.1) on Ubuntu 14.04 LTS following this guide. Everything works - I can get test OSM layer's tiles through proxy. But in order to use this layer on my web page I have to add additional response header:

Access-Control-Allow-Origin: http://my-domain

If found headers option in the documentation but it's meant to

  add additional HTTP headers to all requests to your sources

but not response headers of proxy itself.

So my question is: how can I add additional HTTP headers to my MapProxy's responses?


texts: QGIS How to label individual sides of polygon

texts:
Using native QGIS layer property definitions, I have successfully labelled each polygon in my layer with its area and perimeter.

Never satisfied, I now wish to label each side of each polygon with its length.

It would seem to be a common-enough requirement but I'm not getting good hits.


texts: Displaying outlines for WebGL polygons in Leaflet, source data coming from GeoJSON object

texts:
I have been using the Plugin Glify for Leaflet. It offers a really great and simple way to use web gl to render maps on Leaflet JS. The problem I'm facing is there are custom options for 'vertexShaderSource' and 'fragmentShaderSource', which require complicated shader code. In the long run it's probably fair I learn to write shaders in javascript, C etc., but in the short run this is something I'm very inexperienced in. Is there a way to draw an outline or border to every polygon drawn? I know this has been done in a 3D context and some other languages, however relating that specifically in this context is seeming more difficult.

The relevant code is fairly basic.

L.glify.shapes({
    map: map,
    data: geojson,
    color: function(index, point) {
        //Some code to choose colours
    },
    opacity: 0.8
});

The output is this, which is coloured polygons features, as a layer onto a map.

In short, my aim is to render borders around these polygons using some sort of shader.


texts: Converting US Census' TIGER files into Tableau-readable format with QGIS?

texts:
I'm a journalist using mapping for dataviz in stories. I've found Carto easy to use, but I'd like to use Tableau for a project involving US Census Bureau's school district TIGER files.
It's easy to get the file and convert to KML. Thought I could then convert to CSV and upload into Tableau, but for some reason, KMLCSV Converter isn't working on my computer. (Don't know why -- I have Java Runtime Environment.)
So my next thought was QGIS. I've installed it and tried to convert a few times according to instructions found online, but I end up with files that don't look right. I'll paste in screenshots of two results I got below.


texts: Which OpenLayers.Protocol should I prefer which design patterns provided at?

texts:
At our scenario, there is a server which stores the feature info (coordinates, properties, etc) into a SQL database. Google Maps are used for a layer. The geoJSON data format are used for communication between our server and web mapping client. We want to create/use a WCF (Windows Communication Foundation) service as a communication protocol. Should I implement OpenLayers.Protocol.WFS as a WCF protocol because OpenLayers supports HTTP or WFS as protocol? Or which OpenLayers.Protocol (HTTP or WFS) should I prefer which design patterns provided at?

Thanks a lot for your helps and clarifications,

Yasemin


texts: Creating different label styles for different attribute in layer using QGIS?

texts:
I would like to create labels with different shapes in the same layer (like ellipse and triangle).

Is it possible?

I just read Creating different label styles for different polygon classes using QGIS?, but it's for the style.


texts: Bing OpenLayers shows OK in QGIS and Composer preview, but is misaligned after exporting to PDF

texts:
I have QGIS 2.2 running on a Mac with iOS Maverick.

I have prepared a map with layers with different CRS, including a Bing OpenLayers as base layer. The project has the 'On the fly coordinate correction' option activated, and it works OK, although sometimes things go out of position, but nothing a refresh from a pan doesn't solve.

Paraphrasing part of question #51986:
"Sometimes when changing the scale (scrolling in and out) the [Bing] layer updates (...) wrong in (...) position compared to the polygon layer. But as soon as I "pan" the map everything reconstructs right. Now the problem I have concerns the print composer."

"The result of any kind of export - pdf, or image or svg - of my map (vector layer + [Bing] layer in the background) differs from what is shown in the preview of the print composer: the vector layer and the Google layer positions are mismatching (...)"

Could this be related to On the fly correction?

Any ideas?


texts:
I need to use the MODIS Terra images, in particular the daily data that contains information about Aerosol Optical Depth. What do i have to do is to obtain the value of AOD in some precise point(in my case are ground monitor stations) and then average this values to obtain monthly data at each monitor stations.
Right now i'm using QGIS to open the image as raster layer, the problem is that when i try to open it, there are lot of layers that are opened at the same time. My question is if some of you have experience with MODIS data, and in particular if you know a better way to open and manipulate the MODIS images.


texts:
I'm looking at coordinate pairs defining Australian maritime boundaries. From here.
The coordinate pairs look like this:
AUS-CS-24 14°05¢09.4168²S 156°38¢26.5411²E ITRF2000

AUS-CS-25 14°05¢16.2817²S 156°38¢56.5798²E ITRF2000

My questions is what does the "¢" mean in the context of coordinate data? Also why is there a "²" at the end?


texts: Configuring value for numberReturned for WFS

texts:
The query

GEOSERVER/wfs?request=GetFeature&version=2.0.0&StoredQuery_ID=all&outputFormat=JSON

gives me the following answer

type    "FeatureCollection"
features
0   {…}
1   {…}
2   {…}
3   {…}
4   {…}
totalFeatures   21350
numberMatched   21350
numberReturned  5
timeStamp   "2019-12-14T12:03:06.475Z"
crs
type    "name"
properties
name    "urn:ogc:def:crs:EPSG::25833"

The relevant part is

totalFeatures   21350
numberMatched   21350
numberReturned  5

My StoredQuery looks like this:

<?xml version="1.0" encoding="UTF-8"?><wfs:StoredQueryDescription xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:fes="http://www.opengis.net/fes/2.0" xmlns:wfs="http://www.opengis.net/wfs/2.0" xmlns:gml="http://www.opengis.net/gml/3.2" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:xlink="http://www.w3.org/1999/xlink" id="all">
    <wfs:QueryExpressionText isPrivate="false" language="urn:ogc:def:queryLanguage:OGC-WFS::WFS_QueryExpression" returnFeatureTypes="KRE_ALKIS:sg_flurstueck_f">
        <wfs:Query wfs:typeNames="KRE_ALKIS:sg_flurstueck_f">

         <fes:Filter>
          <fes:Within>
            <fes:ValueReference>shape</fes:ValueReference>
                        <gml:Envelope srsName="urn:ogc:def:crs:EPSG::25833">
                            <gml:lowerCorner>353038.3125 5840211.5</gml:lowerCorner>
                            <gml:upperCorner>378832.09375 5858933.5</gml:upperCorner>
                        </gml:Envelope>
          </fes:Within>
        </fes:Filter>

        </wfs:Query>
    </wfs:QueryExpressionText>
</wfs:StoredQueryDescription>

I would like to have all 21350 features returned and not only 5. Where can I change the setting for valueReturned value?

What have I tried so far?

1.
I added count=100 to the query.

GEOSERVER/wfs?request=GetFeature&version=2.0.0&StoredQuery_ID=all&outputFormat=JSON&count=100

count=100 is ignored by GeoServer. I see the same result but no error.

2.
I tried GEOSERVER/wfs?request=GetFeature&version=2.0.0&outputFormat=JSON&typeNames=sg_flurstueck_f which led to the same result.

This is the information on the GeoServer that should show me 50 features?

Where can I configure the value for numberReturned?


texts: Seeking free solution for intranet geodatabase for point data, connectable to QGIS/MapInfo?

texts:
I have a plenty of point measurements - thousands of points with date/time, gps, measured values etc. I would like to make it possible to put all my particular datasets (can decide whether shp or something like CSV) somewhere and then to be able to load the points in QGIS (my prefered app) and possibly also to newest MapInfo Pro (colleague uses it) on PCs in our local network - but not always all the points but for example particular date range only.

I will have one linux desktop (probably *buntu) which I can use as server and already tried PostGIS some time ago. My questions are following:

is it better to use 32bit or 64bit (multiarch) *buntu? I do not know whether 64bit architecture brings some advantage here...
is PostgreSQL/PostGIS good solution or should I look for something else?

according to the internet - even MapInfo is able to read the data from it

can I make filtered data exports from the database in something like csv or shp?
is it possible to define some script loading a bunch of files (same structure and CRS) into DB just by running it? (looks like I could use this maybe)


texts: Why aren't there LiDAR records in some black and white parts of the scene?

texts:
This is the actual scene:

This is the corresponding LiDAR scene:

The LiDAR data is from xyz file. The questions are:

Why are there no points for the dark squares on the checkerboard as well as the wall behind? Similarly, why are there no points for the floor but points for the "holes"?
Are xyz files preprocessed so that the undesired points are "removed"?

I'm new to LIDAR and point cloud processing so I'm seeking for help here.


texts: Batch conversion of Modis Swath (level 2) data

texts:
I'm using very MODIS Swath products and need to import some bands to TIF. I'm working with MOD06_L2, MOD07_L2 and MOD11_L2 products and the MRT Swath isn't work as well. I'm using the following code:

gdalwarp -of GTIFF -tps -t_srs EPSG:4674 HDF4_EOS:EOS_SWATH:"MOD11_L2.A2014001.1305.006.2016179220648.hdf":MOD_Swath_LST:LST LST.tif

Unfortunately, the output image is blank. I need to convert 3 years of data for the products listed above.

Can anyone help me?


texts: Detecting onClick event for clicking Printout button of ArcGIS API for JavaScript?

texts:
I am wondering if there is a way to detect when a user clicks the "Printout" link after using the print task in ArcGIS Javascript 3.14.  I have set up many custom options in my mapping application where I ended up using a "Prep Map" button to recreate the print widget based on the custom options.  Rather than the Print button resetting to the original state after clicking the Printout link, I would like to hide the print button so need to figure out how to use a function for a printout button click event in order to do so.  Below is the code I tried after using the HTML selector and thinking I could use the class ".esriPrintout" to get the click event but nothing happens.

function openPDF(){
    $('.esriPrintout').on('click', function(){
        console.log("PDF Opened");
    });
}


texts:
It my be old issue, but I can't find right solution.

How to produce GeoTIFF from layout view not from data view?

Any simple solution?


texts: Using cursor.next() is resetting "row in cursor" causing it to skip every other row

texts:
I'm trying to update a column based on the value in the next row, but cursor.next() is giving unexpected behavior. When I identify the next_row using cursor.next(), it is resetting my row in cursor! The result is my loop goes over every other row, not every row. What am I missing? (Code below is just to show cursor behavior; code for row update for column has been stripped.)

Code

import arcpy
temp1 = r'mypath\mypointfile.shp'
with arcpy.da.UpdateCursor(temp1, ['OID@']) as cursor:
    for row in cursor:
        try:
            row_next = cursor.next()
            print 'current_oid=', row[0], 'next_oid=', row_next[0]
        except StopIteration:
            print 'no rows'
        except:
            print 'something else went wrong'

Actual Prints

current_oid= 0 next_oid= 1
current_oid= 2 next_oid= 3
current_oid= 4 next_oid= 5
current_oid= 6 next_oid= 7
current_oid= 8 next_oid= 9

Expected Prints

current_oid= 0 next_oid= 1
current_oid= 1 next_oid= 2
current_oid= 2 next_oid= 3
current_oid= 3 next_oid= 4
current_oid= 4 next_oid= 5

UPDATE

Purpose of this is to only keep rows that are the first or last point in the group.

Example attributes:

oid group   id
0   1   1
1   1   2
2   1   3
3   2   1
4   2   2
5   3   1
6   3   2
7   3   3
8   3   4
9   3   5

Desired attributes:

group   id
1   1
1   3
2   1
2   2
3   1
3   5

Original code, that didn't work because cursor.next() skips rows:

with arcpy.da.UpdateCursor(temp1, ['OID@', 'group', 'id']) as cursor:
    for row in cursor:
        group_current = row[1]
        id_current = row[2]
        if id_current == 1:
            check = "keep"
        else:
            try:
                row_next = cursor.next()
                group_next = row_next[1]
                id_next = row_next[2]
                if id_next == 1:
                    check = "keep"
                else:
                    check = "remove"
            except StopIteration:
                check = "keep"
            except:
                print 'something else went wrong'
            if check != "keep":
                cursor.deleteRow()


texts: Convert coordinates in a CSV file and get the new ones added in the text file

texts:
In QGIS I am trying to convert coordinates from a CSV file. QGIS can use another projection, but when I try to export this layer (and changing the output CRS) the coordinates in the CSV file itself won't change. Is there a way to transform/convert these coordinates and get the new values exported?
Update of details:
I have coordinates for points in a circle around a center point (geocoded from a service using WGS84).

These coordinates are in WGS84. I want these converted to SWEREF 99 TM, since that is what my other software is using. But when I import them, it looks like this:

I have use field calculator with on a CSV file and a SHP file, using this formula:
x(transform($geometry,'EPSG:4326','EPSG:3006'))

The attribute form for the CSV file looks like this:

The formula seems to be working.
The SHP file (exported from the CSV file) attribute form looks like this:

The formula is creating NULL values.
Still, if I use the converted coordinates from the CSV file (export it to a new CSV file)
Most likely there is some basic knowledge regarding coordinate transformation that I do not possess...


texts: Convert Eps to geographic distance using DBSCAN

texts:
I am using the dbscan cluster (package fpc) in R to find clusters on a set of latitudes/longitudes coordinates. I want to find an eps distance that corresponds to a meaningful geographic distance (e.g., a mile or kilometer)
My data looks like this:
  Longitude Latitude
1 -87.53163 41.68640
2 -87.59986 41.67341
3 -87.80099 41.95469
4 -87.82481 41.97409
5 -87.67671 41.68832
6 -87.67751 41.73192

Through trial and error, I got a set of clusters that look plausible. In the attached map.
dbscan.clust1 <- dbscan(points, eps=.025, MinPts=20)

What I'd like to do, is find the eps distance that would correspond to a meaningful measurement (i.e., 1 mile = .025 eps)


texts: Greater Than/Less Data Driven Page Definition Query?

texts:
I'm using data driven pages where each page is representing a series of points for a particular year.  The problem I have is that I want to not only show the points for a particular year, I want to show each point before said year as well.  So as I go through the pages, points are only being added to the map, not removed. I want to do a page definition query so that it would say:
 "year" <= DataDrivenPageYear

The problem is that the Page Definition query only offers "Match" and "Don't Match."  Is there something I can type in the regular query builder to tell it that I want the year less than or equal to the Data Driven Page year?


texts: Conversion from Point 2D to Point 3D with measures

texts: I have a number of PostgreSQL empty tables I am trying to setup for my data. I added the Point Geometry to one but later discovered it is X,Y aware and not X,Y,Z aware. How can I change this? I could only find 3D to 2D on the forms and ST_Force_3D didn't work. Also I am going to need X,Y,Z,M capability on some of my geometries; what do I need to do to impose these values to my table geometries?


texts:
I think my problem is related to this and that question.

I have some big raster layers on my GeoServer (v2.11.2) which span the whole world. They are organized in a Group Layer, where each (sub)continent makes up for one single layer.

On startup, the map is zoomed to the maximum extent and shows everything perfectly fine, but when I pan the map until the dateline comes into sight, the data beyond disappears, meaning the data itself is not wrapped around the world, while the basemap is. How can I fix this?

Map Projection: Pseudo Mercator (EPSG: 3857)
Data Projection: Pseudo Mercator (EPSG: 3857)
Bounding Box set from layer bounds (using CRS bounds does not change the results)


texts: Most accurate GPS hardware for less than 1 mile radius

texts:
I apologize if this is a duplicate but I can't find a good answer or area to know where to find the correct answer.

I am looking for the technology or hardware around wearables or similar technology that can detect movement or at least log location within about a 2 foot area.  I have looked around at it looks like the most common devices use triangulation but aren't more accurate than around 400 feet.

Is there any tech that is more specific than a two foot radius?  I am not trying to waste anyone's time but I can't find a better place to get an answer....


texts: Split features of multi-part polygons according to condition(s) in QGIS

texts:
I have a multi-part polygon with a large number of irregularly-shaped features (>600 features). I would like to split them into smaller sub-features/ segments(?) and assign a unique id to these smaller sub-features. Ideally, I would like to split these features based on the condition that the line/ border between them is straight. However, some of the features hardly have straight lines running through them, so another way to split them could be based on a pre-determined area size.

Is there a function to do this automatically in QGIS? And, what are the steps to do this? I'm using QGIS version 3.4.4-Madeira.

For example, I would like to do something like this and then add a unique id row in the feature table:


texts: Qgis Changing title automatically according to selected sublayer using QGIS composer

texts:
I have a group in QGIS with around 20 sublayers. I would like to automatically update the label in my Layout based on the sublayer I have selected, so i don't have to do this manually.
I have found this code
from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom')
def func(feature, parent):
    root = QgsProject.instance().layerTreeRoot()
    layer = qgis.utils.iface.activeLayer()
    lyr = root.findLayer(QgsProject.instance().mapLayersByName(layer.name())[0].id())
    group = lyr.parent().name()
    if group != '':
        return lyr.parent().name()
    else:
        return 'Group not found'

This somewhat does what I want, but it shows the parent name. And I would like for it to show the name of the selected child. I'm a student and just started learning QGIS.


texts: Using GDAL to extract Corner Coordinates based on a jpg file's GPS EXIF data?

texts:
I have ~100 jpeg images taken by a Drone. I tried to use gdal_merge to merge them and create a big tile GeoTIFF image. gdal_merge spends lots of time, and in the end, it just put each jpeg image precisely on top of the previous one. The GeoTIFF result shows the last jpeg file.

I tried to diagnose the problem and used gdalinfo. I found the jpeg files include full GPS EXIF data, but they do not have proper corner coordinates. That means the corner coordinate of all jpeg images shows like below:

Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 3648.0)
Upper Right ( 4864.0,    0.0)
Lower Right ( 4864.0, 3648.0)
Center      ( 2432.0, 1824.0)

Is there any way to use GDAL to extract Corner Coordinates based on GPS coordination in the jpg’s EXIF file?

I know OpenDroneMap uses the GDAL and perfectly merges Drone images. I just don’t know how it extracts Corner Coordinates from EXIF? Below link  is an example (images are in the images folder):

https://github.com/OpenDroneMap/odm_data_aukerman

Per @Spacedman ask, here is the full extract from gdalinfo. The sample image is here:
https://github.com/OpenDroneMap/odm_data_aukerman/blob/master/images/DSC00229.JPG

OpenDroneMap can perfectly merge the images in the sample meantioned above. However, gdal_merge fails since it cannot figure out Corner Coordinates (my guess)

Driver: JPEG/JPEG JFIF
Files: 0.JPG
Size is 4896, 3672
Coordinate System is `'
Metadata:
  EXIF_BrightnessValue=(7.53672)
  EXIF_ColorSpace=1
  EXIF_ComponentsConfiguration=0x01 0x02 0x03 0x00
  EXIF_CompressedBitsPerPixel=(3)
  EXIF_Contrast=0
  EXIF_CustomRendered=0
  EXIF_DateTime=2016:06:29 10:47:04
  EXIF_DateTimeDigitized=2016:06:29 10:47:04
  EXIF_DateTimeOriginal=2016:06:29 10:47:04
  EXIF_DigitalZoomRatio=(1)
  EXIF_ExifVersion=0230
  EXIF_ExposureBiasValue=(0)
  EXIF_ExposureMode=0
  EXIF_ExposureProgram=2
  EXIF_ExposureTime=(0.004)
  EXIF_FileSource=0x03
  EXIF_Flash=16
  EXIF_FlashpixVersion=0100
  EXIF_FNumber=(3.3)
  EXIF_FocalLength=(4.45)
  EXIF_GPSAltitude=(346.315)
  EXIF_GPSLatitude=(41) (18) (13.7576)
  EXIF_GPSLatitudeRef=N
  EXIF_GPSLongitude=(81) (45) (1.6792)
  EXIF_GPSLongitudeRef=W
  EXIF_ImageDescription=
  EXIF_Interoperability_Index=R98
  EXIF_Interoperability_Version=0x30 0x31 0x30 0x30
  EXIF_ISOSpeedRatings=100
  EXIF_LensSpecification=(4.45) (44.5) (3.3) (5.9)
  EXIF_LightSource=0
  EXIF_Make=SONY
  EXIF_MakerNote=SONY DSC
  EXIF_MaxApertureValue=(3.44531)
  EXIF_MeteringMode=5
  EXIF_Model=DSC-WX220
  EXIF_PixelXDimension=4896
  EXIF_PixelYDimension=3672
  EXIF_RecommendedExposureIndex=100
  EXIF_ResolutionUnit=2
  EXIF_Saturation=0
  EXIF_SceneCaptureType=0
  EXIF_SceneType=0x01
  EXIF_SensitivityType=2
  EXIF_Sharpness=0
  EXIF_Software=eBee
  EXIF_UserComment=
  EXIF_WhiteBalance=0
  EXIF_XResolution=(350)
  EXIF_YCbCrPositioning=2
  EXIF_YResolution=(350)
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 3672.0)
Upper Right ( 4896.0,    0.0)
Lower Right ( 4896.0, 3672.0)
Center      ( 2448.0, 1836.0)
Band 1 Block=4896x1 Type=Byte, ColorInterp=Red
  Overviews: 2448x1836, 1224x918, 612x459, 160x120
  Image Structure Metadata:
    COMPRESSION=JPEG
Band 2 Block=4896x1 Type=Byte, ColorInterp=Green
  Overviews: 2448x1836, 1224x918, 612x459, 160x120
  Image Structure Metadata:
    COMPRESSION=JPEG
Band 3 Block=4896x1 Type=Byte, ColorInterp=Blue
  Overviews: 2448x1836, 1224x918, 612x459, 160x120
  Image Structure Metadata:
    COMPRESSION=JPEG

@RoVo should note that this has to be part of a batch process. I cannot use QGIS or interactive tools.

I am aware of so many commercial tools and products that can do this.

I am asking if I can use GDAL to solve this problem.


texts: modifiedmk Mann-Kendall Test Over NC /Raster stack

texts:
I am working with NetCDF file for analysis of trend. I am using the following codes for analysis of trend..

library(modifiedmk)
library(raster)

    I <- brick("I.nc")
    d <- stack(I)
    tsfun = function(x,na.rm){
      if(all(is.na(x))){return(NaN)}
      return(pwmk(x)[["Sen's Slope"]])
    }

    t <- calc(d, fun=tsfun, na.rm = TRUE)

The above codes are doing good but when I'm replacing [["Sen's Slope"]] with other names it is throwing an error saying "Error in pwmk(x)[["P-value"]] : subscript out of bounds"

The function is stable with only two arguments [["Sen's Slope"]] and [["old. Sen's Slope"]]. I need P-Value, Tau and Z-value. Is there a solution to this?


texts: Web AppBuilder's "Near Me" widget doesn't work with my point layer?

texts:
So I added a zipped shapefile of point data from a file.  This was put on a Web Map and this Web Map is connected to a Web AppBuilder.

I would like to use this "Near Me" widget to locate points closest to an address that is entered. However, the widget says my layers aren't configured properly.   How can I configure my layer or make this work?  Am I missing something in the attribute field?


texts: How to Build Query Table from two OLE DB tables using ArcPy?

texts:
When I run the following code to create a subset of a table, the error says "An invalid SQL statement was used" [Query]. The error shows in line 39. Is there any idea you can figure this out?

import arcpy
# MakeQueryTableOLEDB.py
# Description: Create a query table from two OLE DB tables using a limited set of
# fields and establishing an equal join.

try:
    # Local variables...
    tableList = ["Database Connections/OLE DB Connection.odc/ORADBA.DLY_AVG_DATA"]

    fieldList = [["ORADBA.DLY_AVG_DATA.OBJECTID", 'OBJECTID'],["ORADBA.DLY_AVG_DATA.STATION_ID", 'STATION_ID'],\
                ["ORADBA.DLY_AVG_DATA.SENSOR_NUM", 'SENSOR_NUM'],["ORADBA.DLY_AVG_DATA.DMONTH", 'DMONTH'],\
                [ "ORADBA.DLY_AVG_DATA.DDAY", 'DDAY'], ["ORADBA.DLY_AVG_DATA.NO_OF_VALUES", 'NO_OF_VALUES'],\
                ["ORADBA.DLY_AVG_DATA.DLY_AVG", 'DLY_AVG']]
    whereClause = "ORADBA.DLY_AVG_DATA.STATION_ID = 'ADM'"
    keyField = "ORADBA.DLY_AVG_DATA.OBJECTID"
    lyrName = "Query"
    #QueryTable = "C:/RG/Projects/CriticalDataNetwork/QueryTableTest.dbf"
    # Make Query Table...
    arcpy.MakeQueryTable_management(tableList,  lyrName, "USE_KEY_FIELDS", keyField, fieldList, whereClause)

    # Print the total rows
    print arcpy.GetCount_management(lyrName)

    # Print the fields
    Fields = arcpy.ListFields(lyrName)
    for field in Fields:
        print Field.name

    # Save as a dBASE file
    arcpy.CopyRows_management(lyrName, "C:/RG/Projects/CriticalDataNetwork/QueryTest.dbf")

except Exception, e:
    # If an error occurred, print line number and error message
    import traceback, sys
    tb = sys.exc_info()[2]
    print "Line %i" % tb.tb_lineno
    print e.message


texts: import OSM to GRASS not a valid flag error

texts:
I've downloaded an OpenStreetMap data extract from geofabrik.de. I installed the v.in.osm extension. When I run

v.in.osm -l input=country-latest.osm.pbf

it works well enough to identify the available OGR layers. However, when I try

v.in.osm input=country-latest.osm.pbf output=countrypoly type=boundary table=multipolygons

I get the following:

  ERROR: v.in.ogr: Sorry, <T> is not a valid flag
  ERROR: v.in.ogr: Sorry, <u> is not a valid flag
  ERROR: v.in.ogr failed
  WARNING: No data base element files found

I've tried different layers within the .pbf, and different .pbf files, but I always get the same errors. Can anyone suggest what I'm doing wrong? I'm using GRASS 7.2.0 on Ubuntu 16.04.


texts:
I'm trying to create a simple join script in python from an arcgis model, where the only steps are joining a map to a table to produce a map with the table's properties. The attributes joined together are "CNTYFPID" from the shapefile map and "fips" from the textfile. I have tried editing the script produced by arcgis, but it's not working. Where am I going wrong?

# Import arcpy module
import arcpy
arcpy.env.workspace = "C:/Users/Documents/ArcGIS/Testing123.mdb"
arcpy.env.overwriteOutput = True

# Script arguments
Map = arcpy.GetParameterAsText(0)
if Map == 'tl_2009_us_county_albers_USGS.shp' or not Map:
    Map = "C:/Users/Documents/ArcGIS/Python/tl_2009_us_county_albers_USGS.shp" #provide a default value if unspecified

Crosswalk = arcpy.GetParameterAsText(1)
if Crosswalk == '#' or not Crosswalk:
    Crosswalk = "C:\\Users\\Documents\\ArcGIS\\Testing123.mdb\\cao_fips_crosswalk_for_gis" # provide a default value if unspecified)

# Local variables:
Map = "C:\\Users\\Documents\\ArcGIS\\geography\\tl_2009_us_county\\tl_2009_us_county_albers_USGS.shp"
Final = Map

arcpy.MakeTableView_management(Crosswalk, "fips")

#Process: Join Field
arcpy.AddMessage("Performing Join Field")
arcpy.JoinField_management(Map, "CNTYIDFP", Crosswalk, "fips", "")

The python version I am using is 2.7 and the arcgis version is 10.2. This is what it looks like in modelbuilder

  Error message:  "Traceback (most recent call last):   File
  "C:/Users/Documents/ArcGIS/Yo2.py", line 28, in
      arcpy.MakeTableView_management(Crosswalk, "fips")   File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py",
  line 6306, in MakeTableView
      raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Input Table: Dataset
  C:\Users\Documents\ArcGIS\Testing123.mdb\cao_fips_crosswalk_for_gis
  does not exist or is not supported Failed to execute (MakeTableView)."


texts:
I am trying to export a series of square images, but I can't seem to understand what's driving the dimensions of the final file. Here is my workflow:

region = PointOfInterest.buffer(200).bounds().getInfo()['coordinates']
task = ee.batch.Export.image.toCloudStorage(
  image = NaipImagery,
  description = 'square_export',
  fileNamePrefix = 'square_export',
  bucket = 'my-bucket-name',
  scale = 1,
  region = region)
task.start()

I expect the output to be a 400x400 pixel square. But instead it is a random number of pixels, that changes depending on which PointOfInterest I am using. Here are the details of one such exported image:

im.size = (504, 401)
im.tile = [
 ('tiff_lzw', (0, 0, 256, 256), 758, 'RGBA'),
 ('tiff_lzw', (256, 0, 512, 256), 266517, 'RGBA'),
 ('tiff_lzw', (0, 256, 256, 512), 495006, 'RGBA'),
 ('tiff_lzw', (256, 256, 512, 512), 617484, 'RGBA')]

Frustratingly, You can see that the tiles go all the way to 512px in each direction, while the size is far smaller than that. Why?!? Other PointsOfInterest are (503, 402), or (503, 401), or...

The result is that I can't display the image, without getting the following error:

display(im)
ValueError: tile cannot extend outside image

Can anyone help me figure out how to export a square image, that has the right size and tiles that allow it to be displayed?


texts: OpenLayers add white background to image/png geoserver layer

texts:
i created a new wms layer in openlayers. i get the layer from geoserver as  image/png and transparent=true.

openlayers rotate the raster 30 deg (due to its coordinates - it's fine!) but add white background. the white background start to disappear when zooming in.

how i prevent openlayers from adding white background to raster after rotate? (openlayers setteing? geoserver settings?)


texts: Query (Return) All Points In Query Task and Query Object

texts:
Using ArcGIS JavaScript API 3.21 I need to grab all Points on a Server. The code works for me with Query.where clause at

  query.where = "FeederID  = 'RMT001'";

but as I said I need to get all Points. I tried to run the query without any Query.where statement but not getting back anything. can you please let me know how to get all Points?

    function execute () {
       var query = new Query();
       query.where = "FeederID  = 'RMT001'";
       query.returnGeometry = true;
       query.outFields = [ "CreationUser", "FName" ];
       queryTask.execute(query, showResults);
    }


texts:
OGR (GDAL 3.2.2) appears to incorrectly read the SRS of GeoJSON features and I would like to know if there is anything I can do to prevent this.
I have two GeoJSON files and both represent the same feature.
File 1 includes an explicit CRS property:
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }

File 2 does not include an explicit CRS property, which I believe means CRS84 is implied (RFC-7946).
Both files are opened by OGR's GeoJSON driver, loaded into a layer, and have their SRS inspected in the following manner.
driver = GetDriverByName("GeoJSON")
datasource = driver.Open("...geojson")
layer = datasource.GetLayerByIndex(0)
print(layer.GetSpatialRef().ExportToWkt())

Both files report the following WKT (it ends with AUTHORITY["EPSG","4326"]])
'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'

However, if I explicitly create a CRS84 SRS via OGR (using exactly the same text as File 1's CRS property), it reports a different WKT.
crs84_srs: Final = SpatialReference()
crs84_srs.SetFromUserInput("urn:ogc:def:crs:OGC:1.3:CRS84")
print(crs84_srs.ExportToWkt())

This reports the following (does not end with EPSG reference)
'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]'

Question: Why does OGR treat GeoJSON as EPSG:4326 instead of CRS84, and is there anything I can do to force it to respect the explicit or implicit presence of CRS84 in the input file? I have reviewed the driver options but I see no likely candidates. The difference between EPSG:4326 and CRS84 is subtle but important - for example is affects axis ordering when calling layer.GetExtent().


texts:
I have a QGIS 3 plugin that is in QDialog. Everything is functional but i want to change it to QDockWidget and dock it within QGIS. How can i convert it to QDockWidget?

I am using QGIS 3 with Pyqt5.

I have tried changing the default behvior from QDialog to QDockWidget, but that mess up with the styles.


texts: Error 000445 in iterate Resample for rasters in Arcpy

texts:
I'm trying to resample rasters using a foor loop. Something is going wrong but i don't know what... This is my code:

import arcpy, os, math
from arcpy import env
from os.path import split, join
from arcpy.sa import *

env.workspace = r"C:\resample da kernel\raster 11x11\input\morfom\cluster1"
outFolder = r"C:\resample da kernel\output_30"
rasterList = arcpy.ListDatasets("*","Raster")
for inRaster in rasterList:
    print inRaster # just to view progress
    outName = os.path.join(outFolder,str(inRaster))
    arcpy.Resample_management(inRaster, outName, "30",  "NEAREST")

the error i receive is:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000445: Extension is invalid for the output raster format.
Failed to execute (Resample).

any helps?


texts: Set display options for more than one layer in ArcGIS Pro

texts:
I am trying to set display options for more than one layer in the ArcGIS Pro software.
The task is to activate 3 raster layers at the same time in order to create a new pseudocolor map with bands.
I know that such can be done in the ArcMap via the Layer List in the Image Analysis window,

but I cannot seem to find it in the Pro version.


texts: Is it possible to query a MySQL database using the results from map.getBounds() with Google Maps API v3?

texts:
I have a MySQL database that contains a large number of polygons, one for each district on a map. They are all quite complicated and it's unrealistic to simply load them all at once. I've looked into using MarkerManager for Google Maps API v3 and this seems to be working as intended. However I'm not sure how to best go about querying my database for polygons that are in the viewport.

I'm thinking that it would be best to store the center of each polygon in a separate column as a simple lat,lng pair so that a query could be run on that - and then return the full polygon.

What is the best way to go about this?

SELECT * FROM polygons WHERE [the center of each polygon is within the bounds of the viewport]

Is this even possible? Am I going about this the wrong way?

And additionally, is there a simple script I can run on my data to return the center of each polygon?


texts: When to use ArcGIS Mosaics with pyramids, Overviews, and mosaic-cache?

texts:
I am trying to optimize raster display with ArcGIS 10.3.1 Server images services.

I have a number of images (250) in tif format which I have displayed using a File geodatabase raster mosaic. According to the size of the Geodatabase (35Mb) it seems that the mosaics are referenced (Do mosaics always reference data or can they have embedded data too?)

There are now various optimization tools but I have really no idea which are relevant regarding rendering speed for services.

What are the differences between Pyramids, Overviews and caches when dealing with image data in a mosaic and which (if not all) should I use to optimize arcgis server images services?


texts: One polygon layer for an area with different projects in the same layer?

texts:
I was wondering if it is possible (maybe not ?) to make a structure so that I can have diffrient information (projekts) over the same area/point in the same layer. This would be for an overview of the projects geographically over all areas and to see information about the projects. I don't know if there is any other way to make this.

Let's say there is a harbour with 3 basins in the same area. The area is called ex. Stark harbour and incloud this 3 basins. I could mark out this area with a polygon or a point. Now let's say that this basins has maintenance like repration every half year and that is a new projektnumber everytime. Now there are 3 basins in this area and I would like to keep this polygons/points as one layer named by the area, Stark harbour.
I could of course creat a new layer for every basin and keep track of the information separately, but this would create so many layers when there are many harbours with a lot of basins around the country.

If there is any other way to make something like this could you show me or give me a link?

I will just post an example of what I would like Qgis to have and maybe clear up my bad explanation.

First a polygone over the area and then you would click with identify to get a pop up

In the attribute field you would be able to click on the project number to get to a new attribute with information regarding that project


texts: Get field list of input layer using @alg decorator

texts:
I'd like to write a simple algorithm to sum up a numeric column in an input vector layer. I can't connect the FIELD input to the SOURCE layer. I tried to add 'parent' parameter but it is unknown keyword argument, however it is listed in the algfactory.py.
Here is my code:
from qgis import processing
from qgis.processing import alg
from qgis.core import QgsProject

@alg(name='AlgSumCol', label='Sum of a column', group='examplescripts',
     group_label='Example scripts')
@alg.input(type=alg.SOURCE, name='INPUT', label='Input vector layer')
@alg.input(type=alg.FIELD, name='FIELD', label='Input field', parent='INPUT')
@alg.output(type=alg.NUMBER, name='SUM', label='Sum')

def AlgSumCol(instance, parameters, context, feedback, inputs):
    """ Sum up a numeric column """
    print(type(instance))
    source = instance.parameterAsSource(parameters, 'INPUT', context)
    field_name = instance.parameterAsString(parameters, 'FIELD', context)
    features = source.getFeatures()
    total = 0
    for feat in features:
        total += feat[field_name]
    return {'SUM': total}

I would like to have the names of numeric attributes of the INPUT layer in the FIELD to select, but I get an error:
TypeError: 'parent' is an unknown keyword argument.
If I remove parent keyword argument and change alg.FIELD to arg.STRING, it works but I have to fill the column name manually.
How can I connect the two inputs and limit to numeric attributes? I can do it using the QgsProcessingAlgorithm class, but I would like to use @arg shorter code.


texts: Automatic creation of polygons from Feature Class points in ArcGIS for Desktop?

texts:
I have two point feature classes, let's call them A and B. They are well-spaced and the points are non-overlapping:

My objective is to create a topology using these two sets of points by connecting the minimum number of A points together to form polygons each of which contains one B point. The result would look something like this:

I am familiar with how to use the Topology rules within a Feature Dataset to ensure that the polygons do not overlap and that the A points appear on the edges of the polygons, whereas the B points are properly contained.

I am not familiar with any way to automatically calculate/draw the edges of the polygons themselves in ArcGIS. Is there a tool (or sequence of tools) I could use to automatically group the A feature class points into the vertices of polygons and then produce an output feature class of polygons?

If not, then the solution I'll have to pursue will be to write a toolbox script to process the A and B feature classes and generate a new polygon output class. If there is a "quick" way to do this, without having to write my own Python interacting with arcPy, I'd love to know.


texts: I am getting false parsing and un-necessary quotes for 3 numeric fields when exporting table to ascii (CSV). Items are all DOUBLE. "675,654","5,567,000" is the resulting export format which should be .....,675654",5567000,..... I have checked all ITEM property settings at table level but not yet in catalogue. Will do that next. This is with Arc GIS 10.1, Windows 7 Enterprise P1. Everything points to the fields being well defined as numeric Double INT, no decimals. Anyone have problems with the ASCII table exporter doing bad parsing of numerics?


texts: Missing list of algs in QGIS 2.14 Processing modeler?

texts:
I changed my QGIS-version from 2.8 to 2.14. I am very happy that v.net algs (grass7) are available now as part of the processing plugin. But now there is no list of algs inside the modeler - so no possibility of building models that way. If I switch back to the processing-version of QGIS 2.8 it is all working fine, but v.net algs are missing then. I tried to combine both versions of the processing-plugin already, but nothing is working out.

Does anyone have the same issues, and is there any solution for it?


texts: Join GPS points obtained from .GPX file in Python Console

texts:
Having loaded the .gpx file into a new layer from python console, I'm unable to access the gps points in the attributes table. Essentially, i'm unable to join these points to form a polygon (buildings) or a line segement (roads) depending on the number of gps points in a gpx file. It's possible to achieve the same in using the QGIS GUI, but I need to perform the same from the python console (to be implemented in a plugin).

Any help or suggestions would be much welcomed.

Thanks


texts: Calculating skewness and kurtosis with command line GDAL?

texts:
I have calculated histogram bins with the -hist option in gdalinfo.

Is there an option to calculate skewness and kurtosis, using this or any other GDAL tool?


texts: Download Survey123 photos attachment with Python using their GlobalID?

texts:
I need to download photos from Surve123 which are stored on ArcGIS server. Now I can download all photos from survey at once, but I need to download them separetely according on their GlobalID.

Now I'm using this code:

myPortal = "https://www.arcgis.com"
username = "my_user_name"
featureItemID = "myFeatureItemID"
outputPath = r"c:\temp"

myGIS = GIS(myPortal, username)

myFLItem = myGIS.content.get(featureItemID)
attachLayer = myFLItem.layers[0]
OIDField = attachLayer.properties.objectIdField
myRecord = attachLayer.query(where='1=1', out_fields=OIDField)
for r in myRecord.features:
    myOID = r.get_value(OIDField)
    attachmenets = attachLayer.attachments.get_list(myOID)
    for attachmenet in attachmenets:
        myDownload = attachLayer.attachments.download(oid=myOID,
                                                      attachment_id=attachmenet['id'],
                                                      save_path=outputPath)
        print("Downloaded: " + myDownload[0])

Is there some spot in this code where I can specify GlobalID of each photo?


texts: Show time series data with non-contiguous times using Leaflet TimeDimension plugin

texts:
My data is made up of squares forming a grid where the value of each square changes with time. The times are equally spaced with 1 hour intervals, however, the data is only available for a certain number of hours each day (daylight). The data covers 365 days of a year and the time is stored with each feature (perhaps an inefficient way to store the time data).

I have already managed to show the animation for the first day (linked question). However, when loading data for 2 days, the first day plays well but the player doesn't stop at the end of the data for that particular day and keeps incrementing the time where there is no data.

I have tried to use  TimeDimension timeDimensionOptions: times but with lack of usage examples, I fear that I got the syntax wrong. Here is the JavaScript code.

What would be the way to visualize my data with non-contiguous times using Leaflet TimeDimension?


texts:
Is it possible to use COPY on a file full of WKT POINTs without dropping the explicit SRID on the destination column?

I'm trying to insert several million rows of WGS84 points and am using COPY for better performance over individual INSERTs.

Content of foo.txt:

POINT(0 0)
POINT(0 1)
POINT(0 2)

Table:

CREATE TABLE foo (
    geom GEOMETRY(POINT, 4326)
);

Command:

COPY foo FROM 'foo.txt';

Error message:

ERROR:  Geometry SRID (0) does not match column SRID (4326)
CONTEXT:  COPY foo, line 1, column geom: "POINT(0 0)"

Versions:

PostgreSQL 9.2
PostGIS 2.3


texts: Python GDAL VRT reading using 50% of available memory. How to limit this?

texts:
When reading in areas of a VRT with vrt.ReadAsArray(...) I see RAM usage will keep increasing up to about 50%. This is fine but becomes a problem for when running on a computing cluster with huge amounts of RAM available.
Is there a way to limit how much RAM gdal uses?


texts:
I need to create a C# application takes a Google Map Url like this:

https://www.google.it/maps/dir/40.7191523,-74.0123231/40.7087892,-74.010857/40.7082982,-73.9994353/@40.7123687,-74.0111298,15.33z/data=!4m2!4m1!3e2

and export a GPX file.

Can someone tell me how can I achieve this? I don't need the whole script but just the idea or base API methods to call.


texts: OpenLayers getGetFeatureInfoUrl with dragBox in CQL

texts:
I tried to request Geoserver (2.16) WMS with the function getGetFeatureInfoUrl (From openlayers 6+)

I tried this by passing the dragBox extent as a CQL filter.

      dragBox.on("boxend", () => {
      var extent = dragBox.getGeometry().getExtent();
      mymap()
        .getLayers()
        .getArray()
        .map((lyr) => {
          if (lyr instanceof ImageLayer) {
            //get the view resolution
            var viewResolution = mymap.getView().getResolution();
            // get the view projection
            var projection = mymap.getView().getProjection();
            // get the source
            var source = lyr.getSource();
            // get the dragbox center
            var center = getCenter(extent);
            // create the url
            var url = source.getFeatureInfoUrl(
              center,
              viewResolution,
              projection,
              {
                INFO_FORMAT: "application/json",
                feature_count: 100,
                cql_filter: `BBOX(geom, ${extent.join(',')} )`

              }
            );
            if (url) {
              axios.get(url).then(
                response => {
                  if ('features' in response.data && response.data.features.length > 0)
                    console.log(response.data.features)
                }
              )

            }
          }
        })

        ;
    });

I got a good url as expected :

http://localhost:8080/geoserver/Jeru/wms?&LAYERS=MyWorkspace%3Astreets&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=MyWorkspace%3Astreets&LAYERS=MyWorkspace%3Astreets&INFO_FORMAT=application%2Fjson&feature_count=100&cql_filter=BBOX(geom%2C%20218279.07862854004%2C631674.8131103516%2C222630.00331115723%2C635171.091873169%20)&I=50&J=50&CRS=EPSG%3A2039&STYLES=&WIDTH=101&HEIGHT=101&BBOX=219146.67372894287%2C632115.0852508545%2C221762.4082107544%2C634730.819732666

URL decoded :

http://localhost:8080/geoserver/Jeru/wms?
LAYERS=MyWorkspace:streets
SERVICE=WMS&VERSION=1.3.0
REQUEST=GetFeatureInfo
FORMAT=image/png
TRANSPARENT=true
QUERY_LAYERS=MyWorkspace:streets
LAYERS=MyWorkspace:streets
INFO_FORMAT=application/json
feature_count=100
cql_filter=BBOX(geom, 218279.07862854004,631674.8131103516,222630.00331115723,635171.091873169 )
I=50
J=50
CRS=EPSG:2039
STYLES=
WIDTH=101
HEIGHT=101
BBOX=219146.67372894287,632115.0852508545,221762.4082107544,634730.819732666

So far so good : but it return me empty array of features :

{"type":"FeatureCollection","features":[],"totalFeatures":"unknown","numberReturned":0,"timeStamp":"2020-05-17T11:14:19.801Z","crs":null}

But if i go to geoserver and pass the same BBOX in cql query : i got features as expexted

here is the CQL query passed to Openlayers viewer of geoserver

And here is the result in the viewer (it get features!)

Any Ideas ?


texts: QGIS 2.14.1 Processing Toolbox GRASS 7 .NVIZ Not working

texts:
I have a problem with the command GRASS 7 .nviz in QGIS 2.14.1 installed through QGIS-OSGeo4W-2.14.1-1-Setup-x86_64.exe.md5sum

all GRASS commands work exept .nviz / nviz7

images:

I test with both .nviz from grass .nviz / nviz7 from grass7

not know if I'm doing it wrong, but when I run the command, appears complete/finished but without any result or output and not control windows is displayed.


texts: Why do these points send me to the middle of the ocean or middle of nowhere?

texts:
I've been beating my head against the wall trying to map some points using Leaflet + OSM.

I put in some sample points (near Bellingham, WA, USA).  I've checked the JSON with JSfiddle and GeoJsonlint.  It's valid.  But they both send me to the middle of the ocean or to the bottom of the map.  Any idea what I'm missing here?

GeoJson

{
"type": "FeatureCollection",
"features": [
    {
        "type": "Feature",
        "geometry": {
            "type": "Point",
            "coordinates": [
                48.759553, -122.488225
            ]
        },
        "properties": {
            "name": "ABBOTT NEIGHBORHOOD PARK",
            "address": "1300  SPRUCE ST"
        }
    },
    {
        "type": "Feature",
        "geometry": {
            "type": "Point",
            "coordinates": [
                48.759553, -121.488225
            ]
        },
        "properties": {
            "name": "DOUBLE OAKS CENTER",
            "address": "1326 WOODWARD AV"
        }
    },
    {
        "type": "Feature",
        "geometry": {
            "type": "Point",
            "coordinates": [
                48.759553, -123.488225
            ]
        },
        "properties": {
            "name": "DOUBLE OAKS NEIGHBORHOOD PARK",
            "address": "2605  DOUBLE OAKS RD"
        }
    },
    {
        "type": "Feature",
        "geometry": {
            "type": "Point",
            "coordinates": [
                48.759553, -124.488225
            ]
        },
        "properties": {
            "name": "DOUBLE OAKS POOL",
            "address": "1200 NEWLAND RD"
        }
    }]}


texts: Setting output path to PostgreSQL in processing function of PyQGIS

texts:
I am studying PyQGIS to analyze gis data.
I want to set output path to PostgreSQL in QGIS 3.8 but I have no idea...

According to processing.algorithmHelp('qgis:simplifygeometries'),
only show OUTPUT: ex> 'postgres:...'

it is my code

mycrs=QgsCoordinateReferenceSystem(3857)
uri=QgsDataSourceUri()
uri.setConnection("localhost","5432","TSDLINK","postgres","password")
geometrycol="geom"
tb_1="tsd_level"
uri.setDataSource("public",tb_1,geometrycol)
layer_6=QgsVectorLayer(uri .uri(False),tb_1,"postgres")
layer_6.setCrs(mycrs,True)
mycontext=QgsProcessingContext()
param_simply={'INPUT':layer_6 ,"METHOD": 0,'TOLERANCE': 3000,'OUTPUT':'??????'}
processing.run('qgis:simplifygeometries',param_simply,feedback=None,context=mycontext)

'output':?????????

How to set path?


texts: Why does Geoserver not serve my primary key/ id?

texts:
I have a postgres/postgis table with the columns:

gid [PK] serial NOT NULL,
geom geometry(MultiPolygon,900913),
crop text

My layer in Geoserver however only shows me the feature types "geom" and "crop", the "gid" is not listed.
Problem is, I need the gid to update data of a particular (the selected) polygon.

Is there a way to get the gid? If not, how do you bypass this problem?


texts: Is it possible to send an OpenLayers map to QGIS for printing?

texts:
I have an existing OpenLayers/GeoServer/PostGIS/GeoExt project and I am looking for a way to print an OpenLayers map that contains WFS layers, vector layers dynamically created by users, and WMS layers.

I currently use MapFish and the GeoServer print plugin to generate maps but MapFish cannot handle the WFS or vector layers.

Is there a way to programmatically call QGIS from javascript and send my exiting OpenLayers map object with all its layers to QGIS so users can print the map?


texts: QGIS Has Rendering Error on JP2000 Raster. Works in ArcGIS

texts:
I've got a large .jp2 raster file that was created using 'Mosaic to New Raster' in ArcGIS. I'm trying to open it in QGIS 2.8.6 (I've tried 2.8.14 as well) and am getting a vertical line that won't render properly. The problem occurs at all zoom levels.

Has anyone else had this issue? Is there a quick fix or is this just a bug in the ecw/jp2 raster plugin?

Edit: Here's what gdalinfo shows:

F:\>gdalinfo wor_wgs.jp2
Driver: JP2ECW/ERDAS JPEG2000 (SDK 5.0)
Files: wor_wgs.jp2
       wor_wgs.jp2.aux.xml
Size is 145095, 107689
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-72.016980560378059,42.486752002435146)
Pixel Size = (0.000003234921015,-0.000003234921015)
Metadata:
  COLORSPACE=MULTIBAND
  COMPRESSION_RATE_TARGET=0
  Corder=PCRL
  VERSION=1
Corner Coordinates:
Upper Left  ( -72.0169806,  42.4867520) ( 72d 1' 1.13"W, 42d29'12.31"N)
Lower Left  ( -72.0169806,  42.1383866) ( 72d 1' 1.13"W, 42d 8'18.19"N)
Upper Right ( -71.5476097,  42.4867520) ( 71d32'51.39"W, 42d29'12.31"N)
Lower Right ( -71.5476097,  42.1383866) ( 71d32'51.39"W, 42d 8'18.19"N)
Center      ( -71.7822951,  42.3125693) ( 71d46'56.26"W, 42d18'45.25"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Undefined
  Min=0.000 Max=255.000
  Minimum=0.000, Maximum=255.000, Mean=55.220, StdDev=56.686
  NoData Value=256
  Overviews: 72547x53844, 36273x26922, 18136x13461, 9068x6730, 4534x3365, 2267x1
682, 1133x841, 566x420, 283x210
  Metadata:
    STATISTICS_COVARIANCES=3213.34265916292,3031.77726231577,2726.40295849355,31
29.54838322411
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=55.220242282501
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=56.686353376831
  Image Structure Metadata:
    COMPRESSION=JP2000
Band 2 Block=256x256 Type=Byte, ColorInterp=Undefined
  Min=0.000 Max=255.000
  Minimum=0.000, Maximum=255.000, Mean=55.180, StdDev=54.204
  NoData Value=256
  Overviews: 72547x53844, 36273x26922, 18136x13461, 9068x6730, 4534x3365, 2267x1
682, 1133x841, 566x420, 283x210
  Metadata:
    STATISTICS_COVARIANCES=3031.77726231577,2938.03869526213,2628.71829197938,31
15.81640621639
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=55.180213092815
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=54.203677875787
  Image Structure Metadata:
    COMPRESSION=JP2000
Band 3 Block=256x256 Type=Byte, ColorInterp=Undefined
  Min=0.000 Max=255.000
  Minimum=0.000, Maximum=255.000, Mean=51.535, StdDev=49.223
  NoData Value=256
  Overviews: 72547x53844, 36273x26922, 18136x13461, 9068x6730, 4534x3365, 2267x1
682, 1133x841, 566x420, 283x210
  Metadata:
    STATISTICS_COVARIANCES=2726.40295849355,2628.71829197938,2422.87224936175,27
24.010779609
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=51.535068822992
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=49.222680233422
  Image Structure Metadata:
    COMPRESSION=JP2000
Band 4 Block=256x256 Type=Byte, ColorInterp=Undefined
  Min=0.000 Max=255.000
  Minimum=0.000, Maximum=255.000, Mean=64.895, StdDev=64.197
  NoData Value=256
  Overviews: 72547x53844, 36273x26922, 18136x13461, 9068x6730, 4534x3365, 2267x1
682, 1133x841, 566x420, 283x210
  Metadata:
    STATISTICS_COVARIANCES=3129.54838322411,3115.81640621639,2724.010779609,4121
.28152183373
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=64.894546059985
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=64.197208053261
  Image Structure Metadata:
    COMPRESSION=JP2000


texts: Internal server error status code 500 ArcMap when adding WMS server

texts:
I'm trying to get access to a WMS server with MapGenie. I have an account set up with a password and a list of URLS, but when I try adding the WMS Server I get a connection error message:

Using Internal server error (status code 500).

Could anyone tell me what the problem might be?
Do I have the wrong URL or is there a problem with the account?


texts: Can QGIS open the file type .pitemx or .pkinfo?

texts:
The following file types are associated with ArcMap(.pitemx) and ArcGIS Pro(.pkinfo).
Can QGIS open the file types .pitemx or .pkinfo?
Please see the example image & weblink below for where to download the file types.
ArcGIS map: https://exchange.maps.arcgis.com/home/item.html?id=9c76293962a041e1853d573ae36c1d34


texts: Trying to use “Pit Remove” tool gives out “'escapeAndJoin' is not defined” error

texts:
I have a similar question to this thread:

Trying to use "Difference" tool gives out "'QgsWKBTypes' is not defined" error

According to image [1], I've received this message when I trying to run TauDEM's Pit Remove algorithm in QGIS 2.18.5:

global name 'escapeAndJoin' is not defined. See log for more details.

I've located the file TauDEMUtils.py [2] but I'm not a programmer so I don't know what to do. Until now, there's no solution to this question.

My System: Ubuntu 16.04.2 LTS

[1] http://bit.ly/2nQfaeY
[2] http://bit.ly/2omCaQg


texts: Generate transit route between two bus stations using GTFS data

texts:
I am trying to create a Python application where I can display the best possible route(s) from one bus stop to another bus stop. I am new to GTFS but I have a basic understanding of the working of the GTFS data. I am trying to figure out how I can obtain the optimal routes along with the shapes of the route between two bus stops to plot in my Folium map.
Can someone give me some direction or links to some working tutorials? How can I manually render the route using the shapes file in the GTFS data?


texts: Looking for literature on correcting one DEM with another

texts: I am looking for some scientific papers (any language) according how to deal with the correction of one DEM with a other (such as ASTER with SRTM). Any idea?


texts: PostGIS nearest points with ST_Distance, kNN

texts:
I need to obtain on each element on one table the closest point of another table. The first table contains traffic signs and the second one the Entrance Halls of the town.
The thing is that I can't use ST_ClosestPoint function and I have to use ST_Distance function and get the min(ST_distance) record but I am quite stuck building the query.

CREATE TABLE traffic_signs
(
  id numeric(8,0) ),
  "GEOMETRY" geometry,
  CONSTRAINT traffic_signs_pkey PRIMARY KEY (id),
  CONSTRAINT traffic_signs_id_key UNIQUE (id)
)
WITH (
  OIDS=TRUE
);

CREATE TABLE entrance_halls
(
  id numeric(8,0) ),
  "GEOMETRY" geometry,
  CONSTRAINT entrance_halls_pkey PRIMARY KEY (id),
  CONSTRAINT entrance_halls_id_key UNIQUE (id)
)
WITH (
  OIDS=TRUE
);

I need to obtain the id of the closest entrnce_hall of every traffic_sign.

My query so far:

SELECT senal.id,port.id,ST_Distance(port."GEOMETRY",senal."GEOMETRY")  as dist
    FROM traffic_signs As senal, entrance_halls As port
    ORDER BY senal.id,port.id,ST_Distance(port."GEOMETRY",senal."GEOMETRY")

With this I am getting the distance from every traffic_sign to every entrance_hall. But how can I get only the minimun distance?

Regards,


texts: Can i use conditional statements in cartodb?

texts:
I want to code my Cartodb.js map so that for each marker, if a column is not empty, it embeds a video into the infowindow, and if the column is empty, then it does noting.I saw this example online and tried to copy it, but it didn't work. http://www.w3schools.com/js/js_if_else.asp

for my code, inside the {{content}} area, I said:
if ({{video_embed_link}} IS NOT NULL) {< iframe width="430" height="242" src="{{video_embed_link}}" frameborder="0"> < /iframe>}

any suggestions on what to do instead?


texts: PostGIS Get the POINT that is x meters away from a fixed point on the route (LINESTRING)

texts:
We have a route network stored in the routes table (LINESTRINGs). Basically every sedment (between two vertexes) is inserted as a separate row in this table.

Among that we have a locations table (type POINT) which contains the coordinates.

What I want to achieve is to get the distance between two points over the road. I have the geometries for points and plenty of linestrings for every road.

What I have done so far after searching is the following query (not working):

SELECT
    ST_Length(ST_Line_Substring(line,
          ST_Line_Locate_Point(line, ST_ClosestPoint(line, fpt.the_geom)),
          ST_Line_Locate_Point(line, ST_ClosestPoint(line, tpt.the_geom)))) As length_m,
          ST_AsText(tpt.the_geom) as to_point,
          ST_AsText(fpt.the_geom) as from_point
FROM
    ST_Segmentize((SELECT ST_LineMerge(ST_Collect(the_geom))
       FROM
           (SELECT the_geom FROM routes WHERE day='2015-12-30') AS sq),  0.00001) As line,
           (SELECT * FROM ways_vertices_pgr WHERE id=3017) AS fpt,
           (SELECT * FROM ways_vertices_pgr WHERE id=18207) AS tpt

The routes represent the routes driven by a car so basically some routes can be used multiple times. Can this be a problem anyhow?

What bothers me are 2 things:
- firstly, how to "merge" plenty of those linestrings to get a full linestring on which I can then get the infromation from functions like ST_Line_Locate_Point where you get (and somewhere set the floating values representing the position on the total lenght of the linestring)
- how to build the query then to get point on the line that is exactly x kilometers away from another point (which is the fixed one for every case). I was playing here to get the distances between two sample points on the line, but sadly cannot move further on.

I have really spent quite some time on this issue and cannot move forward. What is interesting in the above query: If I limit the select from routes to 1, it seems to be working, but of course only one part of the routes is added, otherwise I get the: ERROR:  line_locate_point: 1st arg isn't a line error.

The sample data sets are:
Linestrings representing small fraction of routes:
"LINESTRING(-1.3326397 50.9174932,-1.3319842 50.9166939)"
"LINESTRING(-1.3333297 50.9183351,-1.3332187 50.9182141,-1.3330436 50.9179954,-1.3326397 50.9174932)"
"LINESTRING(-1.3338982 50.9186728,-1.333581 50.9185143,-1.3333297 50.9183351)"
"LINESTRING(-1.3341242 50.9187719,-1.3338982 50.9186728)"
"LINESTRING(-1.335291 50.919197,-1.3346072 50.918916,-1.3341242 50.9187719)"

fpt (from point - the starting point) coordinates:
"POINT(-1.3333297 50.9183351)"

The final goal is to get the tpt coordinates for the point that is on the linestrings (on the route) that is x, let's say 100 meters, away from the fpt over the route.


texts: Is it possible to connect QuickMapServices to Apple map satellite imagery?

texts:
Browsing the web I came upon https://satellites.pro/ and noticed that Apple maps has higher satellite resolution imagery than Google or Bing in some areas.  I tried QuickMapServices plugin in QGIS but Apple maps is not listed anymore.
Any idea?  Tried using satellite.pro itself but it didn't work;
e.g.
satellites.pro/Mali_map#14.357793,-2.786880,19
URL: https://satellites.pro/Mali_map#{y},{x},{z}


texts: I want to change the QGIS logo(like in the image) when QGIS starts. Can me help someone for this?


texts:
I am working with an American translation (1922) of the logs of Vitus Bering's voyages in the North Pacific Ocean (1742).
Latitude and longitude are written in the format 52° 02′ N and 1° 57′ E, but the rhumb lines are written in the format S54°21′E.
Are there spaces missing from the rhumbs or is this how they should be written? Why is there a cardinal direction at both the start and end of the rhumb?


texts: Merge two points of two different polygons to get one polygon at the end

texts:
I would like to join points from two different features. I already searched for plugins, but so far but didn't find one working plugin or solution for my task.

So here is a picture of 2 polygons/feautures:

I already marked 2 points and these shall be merged with the other both on the left side. So at the end I would like to have one polygon and the 4 points in the middle for sure shall be 2.

Similar questions:

Join attributes of point and line layer by location (closest distance) in QGIS
how to connect points to create a polygon boundary in QGIS
How to merge lines with slightly different endpoints in QGIS?


texts: Generating 3D surface between 2 linestrings using PostGIS?

texts:
I'm trying to create a surface with a slope between 2 circles, i.e. one large external circle with a Z equal to 300 (m) and another one inside, smaller and  with a Z equal to 200 (m).  I should obtain a polygonal surface like the section of a cone standing upside down.  The Polygon is generated but when a query the min and max Z value of it I obtain the same value.  Here is my code :

SELECT e.id,
    st_Zmin(e.geom) AS max_alti,
    st_Zmin(i.geom) AS min_alti,
    st_Zmin(st_buildarea(st_collect(e.geom, i.geom))),
    st_Zmax(st_buildarea(st_collect(e.geom, i.geom))),
    st_buildarea(st_collect(e.geom, i.geom)))::geometry(PolygonZ,31370) AS geom
   FROM external_circle e
     LEFT JOIN internal_circle i ON e.id = i.id;

id  max_alti    min_alti    st_Zmin st_Zmax geom
1   332         232         332.0   332.0   -
2   346         246         346.0   346.0.  -

It seems that what I obtain is a horizontal donut.


texts: Is it possible to know last update of data behind a WFS service?

texts:
I need to use several WFS service which are not mine and convert them to shapefiles.
I'd like to update the shapefiles whenever the data is updated in the WFS services without downloading the whole data.
Is there a way to know the moment when data changed?


texts:
I want to do a navigation application OSM PBF files.

I'm using the pbf files from the osm website.

I was looking for a way to access the pbf files, I found that it can be accessed using the osmosis library. But I have no idea on how to use it. and there's no documentation.


texts: QGIS 3.2 Quick Map Services Plugin not Showing Basemaps

texts:
For some reason my landsat imagery through QuickMapServices is not showing up. I set the CRS to Albers Equal Area NAD 5070. It does say "Getting Tiles", so I'm not sure if it just takes hours to load?

Any ideas on how to get the imagery to come up?


texts: R {sp} spatial subset: keep so many polygons as number of overlapping points

texts:
I have a set of overlapping spatial points and spatial polygons. I need to subset the polygons by overlapping points, but, if my polygon contains 2 points, I need to keep polygon twice, etc.

To keep each polygon only once, I can simply use sp::over in following statement:

# keep polygons that contains points
outFc<- polygons[!is.na(sp::over(polygons,
                                  sp::geometry(points))),]

How to keep multiple polygons, if multiple points overlaps it?

My example:

library(sp)
library(rgeos)

# Create points and convert them to spatialPointsDataFrame
pts = data.frame(cbind(c(-0.3, -0.2, 0.7, 0.7, 0.1),
                       c(-0.3, -0.2, 1.9, 0.8, 0.4)))

pointsDF = SpatialPointsDataFrame(pts, pts)
pointsDF@data$IDs = paste0("pts", 1:nrow(pointsDF@data))

# Create multiple polygons and merge them into one file
p1 = readWKT("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", id = "stand1")
p2 = readWKT("POLYGON((0 1, 1 1, 1 2, 0 2, 0 1))", id = "stand2")  # 0.5 1,0 2,1 2,0.5 1
p3 = readWKT("POLYGON((0 -2, -1 -2, -1 -1, 0 -1, 0 -2))", id= "stand3")   # 0.5 0.5,0 1.5,1 1.5,0.5 0.5

polys <- rbind(p1,p2,p3, makeUniqueIDs = TRUE)

# windows()
plot(polys)
plot(pointsDF, add = T, col = "red")

# keep polygons that contains points
#-----------------------------------------------------------
# ??? How to change this to keep one polygon for each point, i.e. multiply the polygons???
# ----------------------------------------------------------
outPoly <- polys[!is.na(sp::over(polys,
                                  sp::geometry(pointsDF))),]

The number of selected polygons is 2, but should be 3. It is ok that my polygons would overlap.

I though that I can:

create new feature class for each row in a points
subset the polygons with my single point
if not emply, add this polygon to final outPoly

something like this (not working yet!):

# Subset each polygon once

outPolyLs <- SpatialPolygons(list())

    for (i in 1:nrow(pointsDF)) {
      row <- pointsDF[i,]
      print(row)
      outPolyLs[[i]]<-polys[!is.na(sp::over(polys,
                            sp::geometry(row))),]
      # check if contains any feature, if yes, add it to final output polygons
    }

But I wonder if it is a correct approach, or there is an easier way?


texts: Making tracking data available to Geoserver

texts:
I have tracking data from GPS, which gives latitude, longitude, heading, etc.

Right now it's stored in a non-standard database (Accumulo) and I'm looking into how best to make the data available for Geoserver to be displayed using OpenLayers.

The data is not in any standard format either, it arrives as part of complex JSON files that report more information.

Essentially I have some text based lat/long data that I can put into anything, though I'd prefer to avoid using PostGIS if possible so I don't have to maintain two databases.

What are my options for getting the data into Geoserver?


texts:
I try to rasterize a multiline data set (streets) from a PostGIS database using QGIS 2.18 and SAGA GIS rasterize tool.

I managed to differentiate between street and not street but I'd need to differentiate the raster value for different road classes (1-10). I tried to set the attribute field as the respective field and chose attribute in the output value field.
Unfortunately, output seems to be more as an index from 1 to 2500 than values from 1 to 10.

I also tried GDAL rasterize and GRASS v.to.rast.value. Output of the GRASS GIS is just binary (street/no street)

In the end I'd need a "categorized" raster for diffeerent kinds of road types and nodata values in the other area. It will be one layer of a suitability analysis which I plan to do with raster calculator.

Question:

How can rasterize a multiline vector file and keep the value of a particular attribute?

update 1: SAGA GIS uses the value of ID field instead of "_objectval_id_.
update 2: GDAL tool would work, but I'm not able to set the raster resolution properly.
update 3: obviously SAGA GIS rasterize only works when the attribute is unambiguous assigned. So I'd have to use indices for every road class or create for every class a single raster layer. correct?


texts: How do I geocode street addresses that have information redacted (e.g. 14XX Main St.)?

texts: My city has an open data section listing crimes and accidents amongst other things. They include approximate time, type of crime, etc. The way they formatted the location Looks something like this "42XX street, city, state, zip", how would I geocode an address that doesn't have a a complete street number? If at all possible. I'm looking to incorporate census tract data as well so I need good spatial accuracy.


texts: Unwanted changes in attribute table after processing

texts:
I am currently having issues with the Eliminate Polygon Parts tool (Cartography Toolbox). I have an Advanced License and ArcGIS Desktop 10.6.1. I run the tool on a polygon feature class, in order to delete features < 20 m2. The resulting polygons look fine. However, there is quite a problem with the attribute table of the feature class - it gets messed up completely.
I have 4 fields (int_300 / tiefe_300 / pq_id / p_typ), containing coded information (short integers). Some of the fields and information was added using the arcpy.AddField_management() and arcpy.da.UpdateCursor() tools. The information in the fields need to remain the same for further processing, but they don't. Somehow the "shape_length" and "shape_area" fields get moved to "pq_id" and "p_typ" respectively, "tiefe_300" just gets lost / set to zeros. See screenshots below.
I have found Shifted columns after Eliminate Polygon Part at GeoNet about a similar issue, but it has no solution on offer.
I suspect that there somehow exists a dependency on the order of the fields... but since I add the fields during processing I cannot do much about the original order.
How can I solve this problem?
Before:

After:


texts:
I recently asked a question and accepted the answer, as it was a good guide on accessing the print composer using python. However, I am still looking for a way to print the custom feature form. Does anyone have ideas on how this can be implemented?

I was thinking of a start by using mapCanvas().saveasImage(), but I'm short on ideas as to how to move forward with a screen capture of the feature form.


texts: Transfer feature details from WMS to raster file

texts:
I need to export a WMS Map in a raster (QGIS) as I would like to do a zonal statistic after I need the WMS raster values. I can see them in the Identify results panel. The WMS has single Band color data but after exporting it, it has 4 bands with only the colour values, not the feature values. (see picture)
Is there a possibility in exporting the WMS to a raster and keeping the values? I always tried it with a tiff. Is another format better?
Does anyone know how to export only one band? this might help already.
If anyone knows how to do it in ArcGIS it would help too.
This are the feature details of the WMS map

These are my export settings

these are the feature details of the exported Tiff with 4 bands

The WMS map Bullerkartan 2012, alla källor: (contextualWMSLegend=0&crs=EPSG:4326&dpiMode=7&featureCount=10&format=image/png&layers=mftemp:luftkvalitet-2015-no2ytor&styles&url=http://kartor.miljo.stockholm.se/geoserver/wms?layers%3Dmfraster:bullerkartan-2003-allakallor)


texts: How do I extract data about the curve of a coast line from a satellite image using Google Earth Image?

texts:
I'm trying to use GEE to analyse the shape and curvature of ocean coast lines from satellite images.  What I've been doing so far is: isolating an image from GEE, saving it as a png or jpg, and loading it into MATLAB for edge detection and curve fitting.

I've looked through GEE's docs, but I didn't see any functions that can do what I'm trying to do.  Does anyone know how I can do these task just using GEE?

Thanks.


texts: undefined symbol: OGR_F_SetFieldNull when installing qgis

texts:
I tried to re-install qgis on Ubuntu 16.04 LTS (xenial) using the ubuntugis repository (http://qgis.org/ubuntugis) but I get the following error messages which I cannot resolve:

Setting up qgis-providers (1:2.18.13+24xenial-ubuntugis) ...
/usr/lib/qgis/crssync: symbol lookup error: /usr/lib/libqgis_core.so.2.18.13: undefined symbol: OGR_F_SetFieldNull
dpkg: error processing package qgis-providers (--configure):
 subprocess installed post-installation script returned error exit status 127
dpkg: dependency problems prevent configuration of qgis:
 qgis depends on qgis-providers (= 1:2.18.13+24xenial-ubuntugis); however:
  Package qgis-providers is not configured yet.

dpkg: error processing package qgis (--configure):
 dependency problems - leaving unconfigured
dpkg: dependency problems prevent configuration of qgis-plugin-globe:
 qgis-plugin-globe depends on qgis (= 1:2.18.13+24xenial-ubuntugis);     however:
  Package qgis is not configured yet.

dpkg: error processing package qgis-plugin-globe (--configure):
 dependency problems - leaving unconfigured
dpkg: dependency problems prevent configuration of qgis-plugin-grass:
 qgis-plugin-grass depends on qgis (= 1:2.18.13+24xenial-ubuntugis); however:
  Package qgis is not configured yet.

dpkg: error processing package qgis-plugin-grass (--configure):
   No apport report written because the error message indicates its a followup error from a previous failure.
   No apport report written because the error message indicates its a followup error from a previous failure.
   No apport report written because MaxReports is reached already

Errors were encountered while processing:
 qgis-providers
 qgis
 qgis-plugin-globe
 qgis-plugin-grass
E: Sub-process /usr/bin/dpkg returned an error code (1)

I've purged QGis and GDAL and re-installed everything. Now, It will be installed! However, it won't start because I encounter the same problem mentioned here: gis.stackexchange.com/questions/234010/qgis-fails-to-load/… But the solution (installing scipy 0.17) did not resolve that problem ...


texts: Getting attribute from overlapping feature from other polygon data using ArcGIS Desktop?

texts:
How do I transfer attribute data from one polygon layer (A) to another (B)?

I am using ArcMap 10.2.

I tried using Identity. It doesn't work.

I also tried converting (A) to point so it has the same attribute. And every polygon of (B) has a point that derived from (A). I just don't know how to copy the point attribute to (B).
By the way, (A) and (B) has no join key.


texts: Create symbology for new layer using Python (arcpy.mapping)

texts:
I have some ArcMap addin that creates new layer and add it into table of content. I need to set up a symbology for the newly created layer based on the attribute table (three differnet values > three diferent colours). Is it possible to manage symbology that way using Python (arcpy.mapping module)?

This is part of my script, but it's only add labes:

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

layer_divide = arcpy.mapping.ListLayers(mxd, "layer", data_frame)[0]
if layer_divide.supports("LABELCLASSES"):
    for lblclass in layer_divide.labelClasses:
        lblclass.expression = '"{}" + [FID] +  "{}"'.format("<FNT size = '20'>","</FNT>")
        lblclass.showClassLabels = True

layer_divide.showLabels = True
arcpy.RefreshActiveView()

I can't find any help how to specify a symbology. Is it even possible?


texts: Creating points with attributes from coordinates in ArcPy

texts:
I know how to create points by coordinates with arcpy using the following code:

import arcpy

ptList =[[20.000,43.000],[25.500, 45.085],[26.574, 46.025], [28.131, 48.124]]
pt = arcpy.Point()
ptGeoms = []
for p in ptList:
   pt.X = p[0]
   pt.Y = p[1]
   ptGeoms.append(arcpy.PointGeometry(pt))

arcpy.CopyFeatures_management(ptGeoms, r"C:\Temp\test.shp")

But in my case, every coordinate pair also has an additional text attribute that needs to be written to the shapefile. What would be the approach?


texts: I am trying to make a floor plan for my office, have managed to make all of the rectangular components (e.g. the rooms) however some of the desk layouts are circular and I haven't been able to work out how to make exact and measured circles. They need to have a radius of 1.45m.


texts:
I have an GeoJSON source with the following specifications:

"properties": {
    "sensorid": 500,
    "Locatie": "Prinsestraat 18",
    "9_13_kw1": 513585.0,
    "13_18_kw1": 806317.0,
    "18_22_kw1": 545899.0,
    "22_9_kw1": 477801.0,
    "9_13_kw2": 519198.0,
    "13_18_kw2": 800637.0,
    "18_22_kw2": 550001.0,
    "22_9_kw2": 514971.0,
    "9_13_kw3": 615974.0,
    "13_18_kw3": 916019.0,
    "18_22_kw3": 562191.0,
    "22_9_kw3": 575459.0,
    "9_13_kw4": 615016.0,
    "13_18_kw4": 967845.0,
    "18_22_kw4": 557373.0,
    "22_9_kw4": 559056.0,
    "Totaal": 10097342
  },
  'geometry': {
    'type': 'Point',
    'coordinates': [x, y]
  }

Every value in a field in this JSON with the format "number_number_string" has to be in a different layer. E.g. 22_9_kw3 has to be in it's own layer, just as 18_22_kw4. The result will look look like this:

As you can see in my example, there are multiple points displayed. Each point has it's own layer, but all the data is in 1 (see above) GeoJSON file.

The only way I know how to do this is to create a lot (16) different ol.layer.Vector's and check the properties during styling. Like this:

var meetpuntenSource = new ol.source.Vector ({
  url: 'points.geojson',
  format: new ol.format.GeoJSON()
});

var meetpuntenObject1 = new ol.layer.Vector ({
  source: meetpuntenSource,
  style: styleFunction
});

var meetpuntenObject2 = new ol.layer.Vector ({
  source: meetpuntenSource,
  style: styleFunction
});

ect.
With styleFunction:

var styleFunction = function(feature, resolution){
    var properties = feature.getProperties();
    if(meetpuntenObject2Toggle) {style = blablabla}

I hope you get the idea.
I want to be able to create different layers from the same source, because otherwise I have to load in a lot of the same data. Is there a better way?

EDIT

To illustrate this further, I explain what you see in my JSON:

["9_13_kw1"] = the "name" of the property
[513585.0] = the value of the property, this value is the radius of the displayed circle (divided by a factor ofc.)

The points represent a sensor-point on the map. The value represents the measured value of that sensor. The more measured, the bigger the circle.

To display all these circles I have to create a layer of EACH of these "name"-value combinations.
Every coordinate will have thus 16 different layers: 4x ..kw1, 4x ..kw2, 4x ..kw3 and 4x ..kw4.


texts: Working with Time Slider and Data Driven Pages

texts:
I have parcel polygon data from 1900 to 2000.  Using the Time Slider tool I can animate the parcel map. I would like to parcel attribute information like grantor name as the map is populated over time. I tried data driven pages but could not get the time slider tool and the dynamic text data drven page name to sync.

Can this be done?


texts: I'm trying to spread my ArcMap legend across 3 columns.  Ideally column #2 and #3 should have their symbology header / layer name above them (see picture). ArcMap however is of the opinion that either one or both headers should be placed at the bottom of the column before. In Legend Properties I tried assigning each layer to it's own column. I tried using them as symbology headers or layer names. What am I overlooking?


texts: Export map with legend using openlayers3

texts: Is there any one who has an idea about how to export a thematic map that is accompanied with a legend (picture) using openlayers3? I searched in the openlayers documentation but I get a map without a legend.


texts: Fixing undesirable lines between tiles from SRTM contour line extraction?

texts:
I'm trying to make a South American topography map based on SRTM data on QGIS.  The problem is a visual one. If the topography contours are extracted from individual SRTM tiles, it results a little straight line connecting each curve from both tiles. If the tiles are merged before and then the curves extracted, it doesn't happen. The ideal would be merge all tiles from South America and then extract the contours, but the program or my computer doesn't handle the process. The archive is too large I guess. It seems not a big problem, but in the end the visual impact of these lines turns to be evident in a larger scale. So:

Is there another way of avoiding these connecting lines from each tile?
If not, can I delete then in a faster way than node by node?


texts: QGIS List Unique Values tool only exports to .html format

texts:
Using QGIS 2.16.2 I would like to export the unique values in a shapefile string attribute field to a new database table (which would contain a single field listing the unique values).  The Vector > Analysis Tools > List Unique Values tool correctly selects the unique values, but the only output option is an .html file.  Is there some way to output to a .dbf, or .csv format?

The shapefile is a line format, so the Dissolve tool does not work (it requires a polygon input), as described in Getting list of distinct values from shapefile field using QGIS?


texts: Conditional "On the Fly" Column with PostGIS

texts:
I have two tables I am working with; tl_2015_us_county which consists of county polygons, and tl_2015_48_puma10 which consists of PUMA polygons. I am currently using the below SQL in order to achieve a sort of conditional assignment of the Within field which denotes whether the county is within the PUMA or vice versa. As an extension, I will need to do a conditional assignment of the Urban field which will be set to true if the corresponding county qualifies (probably based off of population density).

How can I achieve this conditional assignment without some sort of crazy number of nested unions and what not? I am not fluent with PostgreSQL/PostGIS.

SELECT
    *
FROM
    (SELECT
        tx_county.name        AS "County Name",
        tx_county.geoid       AS "GEOID",
        puma.pumace10         AS "PUMA",
        'TRUE'::Character(5)  AS "Within",
        'FALSE'::Character(5) AS "Urban"
    FROM
        tl_2015_48_puma10 as puma,
        (select * from tl_2015_us_county where statefp = '48') as tx_county
    WHERE
        st_within(st_centroid(tx_county.geom), puma.geom)
        AND
        tx_county.aland <= puma.aland10

    UNION

    SELECT
        tx_county.name        AS "County Name",
        tx_county.geoid       AS "GEOID",
        puma.pumace10         AS "PUMA",
        'FALSE'::Character(5) AS "Within",
        'FALSE'::Character(5) AS "Urban"
    FROM
        tl_2015_48_puma10 as puma,
        (SELECT * FROM tl_2015_us_county WHERE statefp = '48') AS tx_county
    WHERE
        st_within(st_centroid(puma.geom), tx_county.geom)
        AND
        tx_county.aland > puma.aland10
    ) as AggregateTable
ORDER BY
    "PUMA", "County Name"

As a side note, this code does work several orders of magnitude faster than my Python Shapely/Fiona implementation...


texts: Teradata Store in Geoserver does not yet recognize Teradata Simplify method

texts:
In my previous question i asked why Geoserver WMS on Teradata is noticeably slow. I looked into Geoserver log and it appears that: GeoServer expects Tessellation for spatial index while Tessellation method is deprecated.

Andrea Aime describes that this is because the code for connection to Teradata store is from 2013 and not kept up with current Teradata implementation.

Now i have another curiosity:

the other reason Teradata-backed GeoServer WMS is slow (well at least compared with default-settings PostGIS) is because GeoServer connection to Teradata does not recognize Teradata newly-implemented SimplifyPreserveTopology method (implemented on March 2014 on Teradata DB version 15.00).

Can anyone confirm the above understanding?


texts: How to combine text and integers in a string field?

texts:
I have field with numbers in it (field type: double, number format: numeric) and I need to create a new field and put text in front of the numbers.  I know that I need to add a new string field, but what is the proper expression to put into field calculator?

For example, if my number is 123, I want to change it to be abc123.

I'm using ArcMap 10.1


texts: I have a GeoTIFF that is in "WGS 84" and unts of degrees, and need to get it in units of meters. How to convert using gdal?

texts:
Below is the gdalinfo, I don't really know enough about GDAL to acomplish changing the units, I also don't want to change the number of cells, which seemed to happen when trying out gdalwarp. But another program I am using (BlenderGIS) needs it in meters. Anybody have a recipe?

Size is 2049, 2049
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-80.031969565414485,26.791587548448831)
Pixel Size = (0.000020129647277,-0.000018050547245)
Metadata:
  AREA_OR_POINT=Point
  TIFFTAG_SOFTWARE=MATLAB 9.6, Mapping Toolbox 4.8
Image Structure Metadata:
  COMPRESSION=PACKBITS
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -80.0319696,  26.7915875) ( 80d 1'55.09"W, 26d47'29.72"N)
Lower Left  ( -80.0319696,  26.7546020) ( 80d 1'55.09"W, 26d45'16.57"N)
Upper Right ( -79.9907239,  26.7915875) ( 79d59'26.61"W, 26d47'29.72"N)
Lower Right ( -79.9907239,  26.7546020) ( 79d59'26.61"W, 26d45'16.57"N)
Center      ( -80.0113467,  26.7730948) ( 80d 0'40.85"W, 26d46'23.14"N)
Band 1 Block=208x208 Type=Float64, ColorInterp=Gray


texts: How to write a vector layer as tab delimited txt file using PyQGIS?

texts:
I'm using the following QGIS 2.8 Python code to save results to a txt file

QgsVectorFileWriter.writeAsVectorFormat(inputLayer, "D:/dstest/my_csv", "CP1250", None, "CSV",True)

I would like to be able to write the results to a TAB delimited txt file as opposed to the comma delimited one produced by the statement above. Can this be done through QgsVectorFileWriter.writeAsVectorFormat or is there another way?


texts: Where/how can I find/buy/make building footprints for Kent County, MI (Grand Rapids, Wisconsin, etc)?

texts:
We need building footprints for Kent County, MI along with their associated street addresses.

Primarily we're looking for Grand Rapids and some of Wisconsin, MI, but if someone had access to this data, we'd like Kent County entirely just to be safe.

Do you know anybody/any company that would be able to provide this, or help us find someone that can?

Any ideas?


texts: I'm trying to extract Sea Surface Temp raster (floating integer) that includes both sea temp and land surface temp.  I want the data outside of a clip.  I've clipped the coastline with a polygon in Arc 10.2, however, the issue is with the remaining desired data outside of a polygon.  I've tried extract by mask without success as it will only return values within the polygon and not outside of it.  Any suggestions would be helpful.


texts:
I did some search on internet and found MDT Professional Edition for Magik programming language but it's not free.

What other Integrated Development Environment (IDE) can be used for running Magik codes?

Is there any free IDE for it?


texts: Get Bounding Box of Map Extent in Stand Alone Application

texts:
I basically have the same question as:
How to quickly calculate bounding box coordinates of a selected polygon with QGIS.
However I am building a stand alone application in python and so do not have the iFace class. I can get to the extent class easily enough with  :

box = layer2.extent()
OR box = self.map_canvas.extent()

Which gives me an object that claims it is an extent, but how do I get the coordinates out?