texts: Dissolving line segments in QGIS if they touch and fall within same class?

texts:
I have line segments that have the following 3 attributes (Class, start_point, end_point).
I would like to merge/dissolve line segments that fall within the same class if they intersect (touch).

i.e. if the end_point of one line segment and the start_point of another line segment are the same (similar values).
(NB. I would like to dissolve lines with the same class value.
The dissolve field can also be a new field which will be generated when streams meet the above condition: ie they fall within the same class as well as the end_point of one segment = start_point of another segment)

The circled numbers represent the class (order) of the streams. Each segment has a corresponding class number as shown. The segments circled in red should be dissolved separately. The same should be done for the segments circled in green. Each segment has a start_point and an end_point. For example for the segments in yellow the end_point of one segment and the start_point of the other coincide.

The lower segment circled in red should not be dissolved with the lower segment circled in green because only their end_points coincide.

How can I proceed to do this in QGIS?


texts: Extract house and road using image Google Earth

texts:
where I live does not have gis data, only images from google earth, so I want to take advantage of this image to build gis data,it is local's road, river and house. I am a beginner with this job. I using "Iso cluster Unsupervised Classification" tool in Arcgis. This is results i get.

Some houses and water are the same color.
Are there any other methods to improve the interpretation results?


texts: QGIS action for inserting text in field by hotlink

texts:
I hope to make a simple action in QGIS, where the field in my layer column "Process" can change from "NULL" to "Finished". I use hotlink to click on the polygon in the layer, where I want the action to be.

I tried to do use type Generic and

[%Process%]=substr('Finished',0,10)

but this did not work.

I hope somebody can help me.


texts: Lambert Conformal Conic (One Standard Parallel) - Choosing 'Latitude of Origin/Standard Parallel'

texts:
From The Lambert Conformal Conic Projection With Two Standard Parallels:

  The standard parallels are usually chosen at one-sixth and five-sixths
  of the total length of that portion of the central meridian to be
  represented.

Are there any similar rules for Lambert Conformal Conic with one Standard Parallel for choosing Latitude of Origin/Standard Parallel for a given region?


texts:
A very simple question: Is there a simple tool in QGIS to graphically visualise geometry formatted as text?

Example geometry:

  "POLYGON((571178 6337246,571178 6402217,598061 6402217,598061
  6337246,571178 6337246))"

Result:

  A layer in QGIS showing the polygon from above.

The intention is to do quick and dirty visualisation for debugging purposes - I feel it is easier to relate to a graphical representation than a coordinate list.


texts:
I am developing a WebApp with Flex ArcGIS and I am trying to display an info window when I roll over a graphic that has been added to the graphics layer. I add a graphics layer when the widget loads and add the required listener for mouse over events, this works because when I try to show an Alert message on mouse over it works but trying to show an info window does nothing. This is the callback I used for mouseOver events

private function onMouseOverHandler(event:MouseEvent):void
         {
             var gr:Graphic = Graphic(event.target);

             map.infoWindow.label = "Info";
             map.infoWindow.closeButtonVisible = true;
             map.infoWindow.show(map.toMapFromStage(event.stageX, event.stageY));
         }


texts: Field format when saving joined layer from xls

texts:
If I perform a join of an xls file to a shapefile and then save it as another shapefile, my numeric fields from the xls are displayed with many decimal ceroes.

Here is an example: The upper table is the result of saving the joined table (which is below).

When I check the type of field in the properties window of the original (joined) shape, it is a double - real- 0 length - 0 precision.
In the new saved shapefile the type of field is double - real- 23 length - 15 precision.

How can I get the numbers with no decimals automatically?


texts: Unable to access and download Sentinel-2 images from Python API with the standard code

texts:
products = api.query(footprint.wkt,
                     date = ('20161101', '20161116'),
                     platformname = 'Sentinel-2',
                     processinglevel = 'Level-2A',
                     cloudcoverpercentage = (0,10)
                    )
print (len(products))

Another Error while accessing sentinel data from APi python for the above same code is:
KeyError    Traceback (most recent call last)
/usr/local/lib/python3.6/dist-packages/sentinelsat/sentinel.py
in _load_subquery(self, query, order_by, limit, offset)

    394     json_feed = response.json()["feed"]
--> 395     if json_feed["opensearch:totalResults"] is None:
    396

We are using some unintended behavior of the server that a null is KeyError: 'opensearch:totalResults'
During handling of the above exception, another exception occurred:
SentinelAPIError  Traceback (most recent call last)
3 frames
/usr/local/lib/python3.6/dist-packages/sentinelsat/sentinel.py in _load_subquery(self, query, order_by, limit, offset)

    401           total_results = int(json_feed["opensearch:totalResults"])
    402       except (ValueError, KeyError):
--> 403           raise SentinelAPIError("API response not valid. JSON decoding failed.", response)
    404
    405       products = json_feed.get("entry", [])

SentinelAPIError: HTTP status 200 OK: API response not valid. JSON decoding failed.


texts: Nearest numpy array element whose value is less than the current element

texts: I'm using ArcGIS 10.4x and Python/numpy for this question. Can anyone suggest an efficient way to determine the array location of the nearest element whose value is less than the search element? As an example, if I provide the script a row/col number of a numpy array (the element value for example = 42.5), what is the most efficient way to find the row/col values for the nearest element whose value is less than the 42.5 value? I will have many source cells to iterate through and more often than not, an element with a lower value will be relatively close, so ideally I'd have a way where the search originates from the elements nearest to the source and radiates outwards until a solution is found.


texts: Define Fill Value for NoData areas during/after reprojection

texts:
I'm using rasterio to reproject rasters. The following code snippet should show roughly what I'm trying:

profile = {'crs': proj4_str,
                   'nodata': np.nan,
                   'dtype': str(dem.height.encoding['dtype']),
                   'count': 1,
                   'transform': dst_transform,
                   'interleave': 'band',
                   'driver': 'GTiff',
                   'width': nx,
                   'height': ny,
                   'tiled': False}

        with rasterio.open(dem_out, 'w', **profile) as dest:
            dst_array = np.empty((ny, nx),
                                 dtype=str(dem.height.encoding['dtype']))
            dst_array[:] = np.nan

            reproject(
                # Source parameters
                source=dem_arr,
                src_crs=dem.pyproj_srs,
                src_transform=src_transform,
                # Destination parameters
                destination=dst_array,
                dst_transform=dst_transform,
                dst_crs=proj4_str,
                # Configuration
                resampling=resampling)

            dest.write(dst_array, 1)

Not caring about the values of single parameters, I run into the problem after reprojection there are areas on my grid which are covered with zeroes. However, I would like to have them covered also with NaNs. That's why I'm saying in the profile that nodata should be np.nan, and fill the dst_array also with NaNs. But it doesn't help:

If I sample the dark blue area, it has the values zero. Does anybody know how I can change this?


texts: Problem with spatial join of large shapefiles

texts:
I want to combine the attributes of five shapefiles into one shapefile using QGIS. Therefore, I used several data processing steps including clipping the five shapefiles to the boundaries of my investigation area using Vector > Geoprocessing Tools > Clip and using Fix Geometries to avoid invalid polygons.  The polygon shape and number of the shapefiles are different, but they have the CRS.
I used Vector > Data Management Tools > Join Attributes by location to join the shapefiles.
But after I joined the third shapefile, the shapefile size has already reached around 1 GB. I have read that the maximum shapefile size is 2 GB. While joining the fourth shapefile, the computation time was really long and my computer crashed around half of it. I tried it a second time, but it happened again.
My three shapefiles look like this:

Temperature_Precipation_Slope with 40.423 Features
(already joined three shapefiles)

Field       Type      Length     Precision
Precipitat  Integer   4          0
Temperatur  Integer   2          0
Slope       Integer   1          0

2. BUEK200 with 101 features
Field       Type      Length     Precision
TKLE_NR     Integer   6          0

3. CLC10 with 2675 features
Field       Type      Length     Precision
clc         Integer   3          0

I already decreased and aggregated my data, by using ranges in the attribute tables to minimize the number of polygons within one shapefile. Moreover, I changed the type and size of the fields with Refactor fields.
I could add the fourth shapefile and save it as a GeoPackage, but the size of the GeoPackage was already 17 GB.
While adding the fifth shapefile and saving the result as a GeoPackage, my laptop crashed again. I guess I don't have enough memory space to save such a large file.
I now try to use PostGIS (which I haven't used before) and managed to upload my shapefiles into the two tables Temperature_Precipation_Slope and Buek200_CLC10.
If I understood correctly I can't use ST_Union, as the shapefiles don't have a common attribute.
SELECT * From "Temperature_Precipation_Slope"
SELECT * From "Buek200_CLC10"

However, I don't understand how to perform the spatial join by location


texts: QGIS different results between $length and measuring in the map

texts:
I use QGIS Lisboa and need the length of some shapefiles. When I use the field calculator the results are always around 30% higher than in real, also than a measuring in the map.

Does anybody have a similar problem or knows about a solution?


texts: Creating pseudo-heatmap based on time spent in AOI using ArcGIS Desktop

texts:
I wanted to see if there was a way to create a heatmap based on time spent in a certain AOI.  What I mean by this is a layer you can toggle on/off that would go from clear to red as you spend more time looking at a certain point.

My reason for this is a good number of people at my company struggle to know where they've already been when we're working with large AOIs and large data sets; it can't get tricky to make sure you've reviewed such a large swath of imagery.
Ideally this layer would be clear but as you move around in the imagery it would make darker patches where you've been.  Maybe this could be based on cursor movement; however, I don't think that's the best option.

I don't actually know whether or not this is feasible but I thought I'd ask.


texts:
I have a whole series of GeoTiff tiles, each one about 250mb in size, in total just over 90gb. I want to be able to display these as a wms into QGIS to 20+ GIS users and in the end make them into a tile cache. My problem is that I don't have enough server space to stick the GeoTiff tiles on to in their current state, so they are living on an external hard-drive.  My server contains GeoServer and I've also got access to MapServer on the same server but I'm less fluent its configuration than I am with GeoServer.

Is there a way that I could reduce the size per GeoTiff to make them more manageable or should I convert them into another file tile, png/jpeg2000 for example?

Has anyone else any experience of such large GeoTiff sizes and what did you do with them?

Oh and I'd like to do this entirely with OpenSource (if possible).

Med


texts: Qgis 3 - Rule-base labeling combining to backgrounds

texts:
I have a pointlayer with no symbology, but is labeled by numbers 1-74. The labels is styled with a circle background. For some of the points I need to keep the circle, but add a small SVG-symbol to the lower right of the label.

I have no idea have to solve this, because as far as I see, I can only choose the circle OR the SVG file.

I have tried to make more labelrules to combine, but it only shows one of them.

The first rule on the pic below "Qr-kode"  =  'Nej' is the standard, and is working fine. The next rule  "Qr-kode"  =  'Ja' and "text = nr" is identic to the first, but goes for the points for which i need to add the SVG file. The last one where "text = tom" is the SVG file. But it only shows either the cicle label or the svg depending on zoomlevel.

Can anyone please help me. Thanks.


texts: OpenLayers map initialised but not displaying

texts:
I am using OpenLayers with Google maps. I have created my map object successfully - there are no errors raised and console.log(map) shows that it is indeed an OpenLayers class - but nothing is showing on screen. All I am seeing is a white rectangle.

Does anyone know what might cause this?

Any tips for finding the bug that is causing this?

With no errors or obvious problems I don't know where to begin. It should be working!


texts: Convert data {shape+XML} format to my custom map data format {XML}

texts:
I work with "map" represented logically as custom XML format not compatible with any OGC format like GML, but roughly based on OSM XML. Unfortunately I cannot change to work with OGC compatible known XML format.

I have data from GIS supplier: shapefile & corresponding XML, also not compatible with any OGC format (of course the shapefile is valid known format).

My task is to convert supplier data to my XML format. It's understood I have attributes not existing in the supplier's data and vice versa.

Because direct conversion from first XML format to second XML format (both are custom not well known formats) means complex piece of code basically merging XML schemas (am I correct here?). I think to apply 2 step approach - First convert to OSM and then from OSM to my format - Which again will result coding piece of code transforming known XML format (OSM) to my XML format, which will be complex task, but less complex that direct conversion (I can take advantage over exiting OSM open sourced codes). (again - am I correct here - it's inevitable down the road code some XML merger algorithm?)

So If I just try to export the shapefile to OSM (with lot's of tools out there) - It's not 100% guarantee all the data from the corresponding XML of the shapefile will be converted also, am I correct? It depends whether the supplier while generating the shapefile included into them all the data from the XML...


texts:
I'm trying to prevent a user in edit mode from creating overlapping shapes. With a CREATED event, I can easily remove the layer. But I want to revert an edited shape back to its original position instead of removing the layer.

My code is as follows:

map.on(L.Draw.Event.EDITED, function (e) {
   var layers = e.layers;
   layers.eachLayer(function (layer) {

     var reversedCoordinates = [], latlngs = layer.getLatLngs();
     for (var i = 0; i < latlngs.length; i++)
         reversedCoordinates.push([latlngs[i].lng, latlngs[i].lat]);

     reversedCoordinates.push([reversedCoordinates[0][0], reversedCoordinates[0][1]]);
     if (checkRectangleOverlap(reversedCoordinates)) {

         //drawnItems.removeLayer(layer);
        //revert the layer to its original position here
      };
    });
 });

Could anybody help me in this regard?


texts: How to randomly select n points within a radius in QGIS?

texts:
I have a dataset of 5000 records which I've plotted over an OSM layer using QGIS and I wish to randomly select 50 from within a radius of 2km from a point of my choosing.

It is similar I think to this here where again the Select Features by Radius is almost what I'm aiming for, but unlike the previous question, I only wish to randomly select 50 points from within the 2km radius.

Using the Select Features by Radius, it looks like there would be approximately 300-400 points within any 2km radius. So a solution would have to be able to randomly select only a fraction of these.


texts: Pansharpening changes background value from 0 to a different value, but only sometimes?

texts:
I have a python script that takes a folder full of multispectral worldview tifs and pansharpens them with their matching panchromatic band.

Both these datasets have a background value of 0 but after I pansharpen it certain dates will have a background value of 1 (sometimes a higher value like 5), which I do not want. I need the value to stay zero for other processes later on.

Basically I am using the Gram-Schmidt method to pansharpen R, G B, NIR2 and then pansharpen the other four bands (RedEdge, Yellow, Coastal blue, NIR1) and then creating a composite to combine the two  pansharpened outputs back into one raster. I do not want to use a mosaic dataset as these outputs are used later outside of ArcGIS.

Here is the code I use:

def pansharpening(multispec_list, pan_list):

    # Pan sharpen each image individually
    if len(multispec_list) != len(pan_list):
        sys.exit("UNEVEN NUMBER OF PAN AND MULTI BANDS")

    pansharp_list = []
    for i in multispec_list:
        # find the unique ID in the multispectral tif and match it to the pan band

        multi = i

        pan = multi.replace("M2AS", "P2AS")

        if pan not in pan_list:
            sys.exit("Missing matching pan bands to " + i)

        # if the ID does not exist in the pan list the variable will be empty after that loop

        pansharp_rgbn2 = os.path.splitext(multi)[0] + "_pansharpen_rgbn2.tif"
        pansharp_recbyn1 = os.path.splitext(multi)[0] + "_pansharpen_recbyn1.tif"
        comp_out = os.path.splitext(multi)[0] + "_pansharpen_allbands_0.tif"
        comp_out_nodata = os.path.splitext(multi)[0] + "_pansharpen_allbands.tif"

        # pansharpen R, G, B, NIR2 and CoastalB, RE, Y, NIR1 bands

        print("pansharpening " + multi + " with " + pan + "and WorldView-3 inputs" )
        arcpy.CreatePansharpenedRasterDataset_management(multi, 5, 3, 2, 8, pansharp_rgbn2, pan, "Gram-Schmidt", "", "", "", "", "WorldView-3")
        print("pansharpened R,G,B,NIR2")
        arcpy.CreatePansharpenedRasterDataset_management(multi, 6, 4, 1, 7, pansharp_recbyn1, pan, "Gram-Schmidt", "", "", "", "", "WorldView-3")
        print("pansharpened CB,Y,RE,NIR1")

        # create a composite with all 8 bands in the correct WV3 order
        coastal = pansharp_recbyn1 + "\Band_3"
        blue = pansharp_rgbn2 + "\Band_3"
        green = pansharp_rgbn2 + "\Band_2"
        yellow = pansharp_recbyn1 + "\Band_2"
        red = pansharp_rgbn2 + "\Band_1"
        rededge = pansharp_recbyn1 + "\Band_1"
        nir1 = pansharp_recbyn1 + "\Band_4"
        nir2 = pansharp_rgbn2 + "\Band_4"

        print("creating composite")
        arcpy.CompositeBands_management ([coastal,blue,green,yellow,red,rededge,nir1,nir2], comp_out)

        # set nodata to 0
        print("setting nodata to 0")
        arcpy.CopyRaster_management(comp_out, comp_out_nodata, "", "", "0", "NONE", "NONE", "16_BIT_UNSIGNED", "NONE", "NONE")
        pansharp_list.append(comp_out_nodata)

        arcpy.Delete_management(pansharp_rgbn2)
        arcpy.Delete_management(pansharp_recbyn1)
        arcpy.Delete_management(comp_out)

I can take the raster and reclassify the 0's to 1's later, but I want to avoid this from happening in the first place. As I'm not sure how to check in python if only the background value is 1s or 0s besides bringing it into ArcGIS Pro and visually looking at the background. (checking the whole raster for 1's in python may not work as I would expect there to be a value of 1 somewhere in the tif).

Any idea why this is happening?


texts: How to find reasonable locations (location study)?

texts:
My GIS skills got a little rusty and I am new in using QGIS (I have learned on ARC).

I have created a map with population estimates and points. I created a buffer around the points of 75km and calculated the population within the buffers (There are overlaps that are not elegant, as they are multiple-counted but thats not my priority). I basically followed this introduction: Calculate number of people within 50 km of a point

Is there a way to calculate potential areas (in this case within europe), where my buffers miss out on largely populated areas?

Basically I would like to receive 75 km buffers located at those areas.


texts: Has anybody successfully added NFHL (National Flood Hazard Layer) to ArcGIS for flex?

texts:
I keep getting the following error when I add the map service to my application:

apServer layer failed to load: Fault code: Channel.Security.Error

Fault info: Security error accessing url

Fault details: Destination: DefaultHTTP

Any ideas?


texts: Leaflet Tooltip bug when added leaflet-draw

texts:
Here's a simple demonstration:
The file leaflet.draw-src.js was downloaded here: https://github.com/Leaflet/Leaflet.draw/tree/master/dist

If you remove the script tag with the leaflet.draw-src.js the tooltip work just fine, but when added the leaflet.draw it gives the error Cannot read property 'popupPane' of undefined.

ps: Not using leaflet-draw is not an option

<!DOCTYPE html>

    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.0.0-rc.3/dist/leaflet.css" />


texts: Python-gdal create geotiff from array with colormapping

texts:
I am creating grayscale, single channel geotiff files from numpy arrays with the below module, but I would like to apply a fairly standard colormap to the data.  I've tried three channels of R,B and G and have had no success if producing anything that looks correct.  I'm thinking my approach, in general is flawed.

How can I go from grayscale to color geotiffs?

def geotiff_output(outfile, mag_grid, lons, lats):
    #*********************Test***************************
    from osgeo import gdal, osr, ogr
    xres = lons[1] - lons[0]
    yres = lats[1] - lats[0]
    ysize = len(lats)
    xsize = len(lons)
    ulx = lons[0] - (xres / 2.)
    uly = lats[-1] - (yres / 2.)
    driver = gdal.GetDriverByName('GTiff')
    ds = driver.Create(sc_settings.dynamic_folder_script + outfile[:-4] + ".tif",
                       xsize, ysize, 1, gdal.GDT_Byte, )

    # this assumes the projection is Geographic lat/lon WGS 84
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    ds.SetProjection(srs.ExportToWkt())
    gt = [ulx, xres, 0, uly, 0, yres ]
    ds.SetGeoTransform(gt)
    outband=ds.GetRasterBand(1)
    #mask no data points...
    mag_grid1 = np.ma.masked_where(mag_grid==-99.,mag_grid)
    #assign actual data values for inclusion in metadata
    outband.SetStatistics(np.round(np.min(mag_grid1),decimals=4),
                            np.round(np.max(mag_grid1),decimals=4),
                            np.round(np.mean(mag_grid1),decimals=4),
                            np.round(np.std(mag_grid1),decimals=4))

    #normalize the array to a 0-255 scale
    if np.min(mag_grid1) < 0:
        mag_grid1 = mag_grid1 + np.abs(np.min(mag_grid1))
    #Normalize array to a 0-255 scale for raster channel
    mag_grid2 = ((mag_grid1 - np.min(mag_grid1))/(np.max(mag_grid1) - np.min(mag_grid1)))*256
    outband.WriteArray(mag_grid2)
    ds = None


texts:
I was wondering if I could use my android phone (HTC Wildfire) to plot gps points and export these into GIS software (I have QGIS and MapInfo).

Thanks,
Alan


texts: .mdb file and MinGX, MinGY, MaxGx, MaxGY parameters

texts: I have .mdb file ( Microsoft Access file ) which appears to be a shapefile index ( snapshot as attached ) and i need to get longitude and latitude from this so that file can be fed to Tableau or PowerBI for visualization purpose.


texts: Basin shapefile into HEC-HMS coordinate system issue

texts:
I am trying to import a basin shapefile into HEC-HMS without going through with the whole HEC-GeoHMS process in an attempt to use HEC-HMS to calculate PMP values for the drainage area.

I make a new HEC-HMS file.
I set the units to US Customary.
I create a basin model.
I define the coordinate system for the basin model using the text provided in the ESRI .prj file.
I import the shapefile of my drainage area (same shapefile I pulled the .prj file text from.
I get the error "The coordinate system of the GIS file does not match the coordinate system of the basin model".

The goal is to use Vermont State Plane NAD83 Feet (2011) as the coordinate system. The shapefile was created / projected in ArcMap.

Below is the .prj file text I am trying to use. I was under the impression that this was the equivalent of WKT.

PROJCS["NAD_1983_2011_StatePlane_Vermont_FIPS_4400_Ft_US",GEOGCS["GCS_NAD_1983_2011",DATUM["D_NAD_1983_2011",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",1640416.666666667],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-72.5],PARAMETER["Scale_Factor",0.9999642857142857],PARAMETER["Latitude_Of_Origin",42.5],UNIT["Foot_US",0.3048006096012192]]

I have spent a fair bit of time trying to solve this issue trying other coordinate systems etc.. I am new to HEC-HMS and use GIS software for very specific purposes; by no means am I an expert user.


texts: Line_interpolate_point: 2nd arg isnt within [0,1] in PostGIS

texts:
Im trying to interpolate a point on edge using ST_LineInterpolatePoint but it is giving an error saying that Line_interpolate_point: 2nd arg isnt within [0,1]'

create table interpolate as (
select coalesce(
        case
        when abb.node=edg1.start_id
        then ST_LineInterpolatePoint(st_makeline(st_linemerge(edg1.geom)), (300-abb.cost)/edg1.travel_spped :: double precision)
        else ST_LineInterpolatePoint(st_makeline(st_linemerge(edg1.geom)),((300-abb.cost)/edg1.travel_spped) :: double precision)

        END,
        abb.geom
        )as geom from  abb left join hyd_network as edg1
        on abb.node in (edg1.start_id,edg1.end_id)
AND edg1.gid not in (select edge from abb) group by abb.node,edg1.start_id,abb.cost,edg1.travel_spped,abb.geom  )

But if i replace 2n argument with any value between 0,1 like ST_LineInterpolatePoint(st_makeline(st_linemerge(edg1.geom)), 0.3) it is giving result as expected. In my code the 2nd  argument also giving value between 0,1 but it is not taking.


texts:
I have a line dataset that contains some lines where the "end points" are actually found in the middle of the line. In the image below, the nodes are shown in green and one of the end points in red. The line doubles back on itself.

How would I ensure the end points are actually at the end of the line, and remove the overlap?

These lines are both valid, and complex - as they have overlaps, as confirmed using Shapely.

>>> from shapely.wkt import loads
>>> wkt = "LINESTRING (0 0, 2 0, 1 0)"
>>> l = loads(wkt)
>>> l.is_valid
True
>>> l.is_simple
False

Ordering the coordinates would work in some cases:

>>> sorted(l.coords)
[(0.0, 0.0), (1.0, 0.0), (2.0, 0.0)]

But not in others as it would change the geometry of the line:

>>> wkt = "LINESTRING (0 0, 2 0, 1 -1)"
>>> l = loads(wkt)
>>> sorted(l.coords)
[(0.0, 0.0), (1.0, -1.0), (2.0, 0.0)]

Attempts at buffering the line with a 0 buffer, and merging lines have not changed the coordinate order.


texts: WMTS Error when adding Mapbox Layer to QGIS: Failed to download capabilities?

texts:
I've tried this a few different ways and can't quite get it to work using QGIS's add WMTS layer tool. The URL loads fine in a browser but won't download as a tile in QGIS. I don't have that good an understanding of how this works/ permissions stuff either - I'm not even sure what it means by capabilities. Any advice?

Full error text:

Failed to download capabilities:
Download of capabilities failed: Error transferring (URL)
-server replied: Not Found

Mapbox URL:https://api.mapbox.com/styles/v1/cdickinson11/cjkwiea7n0raq2rp9gvab9hc4/tiles/256/{level}/{col}/{row}@2x?access_token=pk.eyJ1IjoiY2RpY2tpbnNvbjExIiwiYSI6ImNqY3J6eDlscTBweHMycXMwcWltbnRiN3MifQ.X68p3M-5wwH3AudbxEYdfw

EDIT: Went back and tried to create a link using this guide from mapbox. Still no luck. New error: not a wms server.

New, homemade link:

https://api.mapbox.com/styles/v1/cdickinson11/cjkwiea7n0raq2rp9gvab9hc4.html?fresh=true&title=true&access_token=pk.eyJ1IjoiY2RpY2tpbnNvbjExIiwiYSI6ImNqY3J6eDlscTBweHMycXMwcWltbnRiN3MifQ.X68p3M-5wwH3AudbxEYdfw#0.8/-0.000000/0.000000/0/wmts?access_token=pk.eyJ1IjoiY2RpY2tpbnNvbjExIiwiYSI6ImNqZThwZGg1dDAyMngzM3Fsdmp0dm82YTcifQ.bjVqusaDZqdZVWunQ0xS7A4


texts: QGIS doesn't handle PostGis connection resets?

texts:
I'm highly frustrated by the fact that if I leave the database alone for a more than a few minutes, I can no longer insert new lines. If I restart QGIS, everything goes back to normal. I believe this is caused by the tcp connection being closed (same problem that causes crashes in PgAdmin)?

Any way to mend that from QGIS side? Or do I have to increase timeouts somehow?

Will pgbouncer help here?

2016-10-05T18:15:38 1   NOTICE: WARNING:  there is no transaction in progress
2016-10-05T18:37:39 1   Query: BEGIN returned 7 [no result buffer]
2016-10-05T18:37:39 1   resetting bad connection.
2016-10-05T18:37:39 1   retry after reset succeeded.
2016-10-05T18:37:39 1   Query failed: SELECT nextval('features_tal_id_seq'::regclass)
            Error: no result buffer


texts: Statistically comparing two sets of coordinates

texts:
I have been given a task to compare two sets of points, one set being a set of predetermined points, on the other being a set of actual points that were taken in the field.

I have to compare these two sets and come up with the statistical difference between the two sets using qgis.

Could someone recommend a process I could go through to achieve this?


texts: ArcGIS JS Search Widget - searching by latitude, longitude

texts:
I have setup a simple map with a default search widget, similar to this:

require([
  "esri/map", "esri/dijit/Search", ...
], function(Map, Search, ...) {
  var map = new Map(...);
  var s = new Search({
    map: map
  },"search");
});

In my web browser, using the widget, I want to search for a location using latitude longitude coordinates (lat,lon), but it seems that the widget requires longitude first (lon,lat) when searching by coordinates. Is there any way to configure the widget to take the lat,lon format? Or would this have to be configured on the server side service that the widget is calling?


texts: Is there any way to match images against google maps satellite view?

texts:
I've got an image that was snapped from google maps satellite view. I know the general location (within a county), but need to find the specific location of the top-down photo.

Are there any tools that can search the google maps satellite images in this given area for a match on the image I have?


texts: Big data and spatial regression: how to make spml more efficient in R?

texts:
I am running a spatial panel regression in R using spml from package splm. When I conduct this on about 1/200th of my data, the program executes. But when I try to do so on a third of my data, R uses up all the memory on the terminal I am working on and the program crashes. My code for specifying the spatial regression is below. How would I make this function operate efficiently enough to produce results and not crash? My terminal has 22 core, so parallel processing is likely an option - I just don't know how to perform this with spml.
results = spml(fm,
                data=dtaSub,
                listw=weightsList,
                model="within",
                effect = "twoway",
               spatial.error = "b")
```


texts: Viewing a small pointcloud in QGIS 3.0 3D Viewer

texts:
I have a pointcloud in postgis, imported using pdal pipeline and chipper to 600 points per patch. I have the pointcloud (only 5000 points) displayed in QGIS. How do I get QGIS to render in in 3D in the 3D viewer, for some reason it is ignoring the z-dimension. It is displayed using the database manager and the following query.

SELECT PC_EXPLODE(pa)::geometry
                 as geom
                 from pcpatches

The screen shot shows the 3D viewer on the right, it should be a deep trench.


texts:
I am a total GIS noob and am trying to read weather data from a PostGIS database. The tables contain data points consisting of a geometry, a timestamp and the actual observed variable. The tables are indexed on geometry and timestamp.

I want to find the values for the observed variables for a specific time range for the nearest point to some coordinates. My grid is in SRID 4283 and a sample query would look something like this:

SELECT * FROM rain
WHERE timestamp >= '2016-11-01' AND timestamp <= '2016-11-02'
ORDER BY rain.geom <-> 'SRID=4283;POINT(141.4539 -31.9539)',
         rain.timestamp OFFSET 0 LIMIT 100

Sometimes this query runs in about a second at other times, it takes minutes. I have not been able to find a correlation with the server load or anything like that.

My tables are reasonably large (2-20 million rows) and I understand that the sort is probably expensive. However, I suspect that some temporary index is built when a query is run because subsequent similar queries seem to be much faster.

Anyway, what I would like to know is how I can make a query with a simple point distance sort like this faster. I know the upper bound of the resolution of my grid. Could I prefilter for a specific boundary box? Would that help with the speed?

EDIT:
The EXPLAIN ANALYZE for the query above looks as follows:

Limit  (cost=239664.44..239664.69 rows=100 width=72) (actual time=175789.131..175789.187 rows=100 loops=1)
  ->  Sort  (cost=239664.44..246555.17 rows=2756291 width=72) (actual time=175789.129..175789.146 rows=100 loops=1)
        Sort Key: ((geom <-> '0101000020BB100000F1F44A5986AE614087A757CA32F43FC0'::geometry)), "timestamp"
        Sort Method: top-N heapsort  Memory: 39kB
        ->  Index Scan using ix_rain_timestamp on rain  (cost=0.43..134320.98 rows=2756291 width=72) (actual time=2.595..174376.437 rows=2710710 loops=1)
              Index Cond: (("timestamp" >= '2016-11-01 00:00:00+00'::timestamp with time zone) AND ("timestamp" <= '2016-11-02 00:00:00+00'::timestamp with time zone))
Planning time: 10.577 ms
Execution time: 175789.282 ms
(8 rows)

This was run on cold caches.

The closest match comes back with a distance of 0.018410813371485 while the furthest match in the limited set is 0.10473267419483. I'm really only interested in the closest point but cannot easily predict how many hits I will get for that one point.

There are 2,710,710 matches for the query without a limit for the timestamps in the WHERE clause. The total number of rows in the table is probably about ~15 million.


texts: ESRI GeoEvent Streaming Services on AWS failing to stream

texts:
I'm attempting to output streaming services on my ESRI Server install on Amazon Cloud but I am unable to.

My Service is successfully being published to server but when I attempt to view it using the default Java Script the stream service URL is red with the error Cannot Connect I opened the chrome development tool and found the error failed: WebSocket opening handshake was canceled

It seems like there is a port conflict between AWS and ESRI but I can't track it down.  I do not have a firewall and I am not restricting any ports. And I know It is not a socket issue because I have bytes set to 1024000 in GeoEvent.


texts: ArcGIS Javascript API Query & Add features without map using REST API

texts:
I need to query arcGIS map service and add this result to another feature service using ArcGIS JavaScript API. I should not use MAP object.

I did query data without map from ArcGIS sample: Now i need to post this data to feature service:

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">

    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
    <title>Query State Info without Map</title>

    US state name :

Any Idea?


texts: How to make the Overlays in OpenLayers exclusives using events at checkbox level?

texts:
On my map I have 3 overlay layers, all of which need to be mutually exclusives, i.e. can't be on (checked) at the same time. I thought of using events to trigger the overlays as the following in my index.html main file:

function init(){
 var map, kml_1Layer, kml_2Layer, kml_3Layer;
...
kml_1Layer.events.register("visibilitychanged", this, layerChanged);
kml_2Layer.events.register("visibilitychanged", this, layerChanged);
kml_3Layer.events.register("visibilitychanged", this, layerChanged);
...
}

function layerChanged(layer) {
     if (layer.object.visibility) {
     switch (layer.object.name) {
        case "kml_1Layer":
              kml_2Layer.setVisibility(false);
              kml_3Layer.setVisibility(false);
              break;
        case "kml_2Layer":
              kml_1Layer.setVisibility(false);
              kml_3Layer.setVisibility(false);
              break;
        case "kml_3Layer":
              kml_1Layer.setVisibility(false);
              kml_2Layer.setVisibility(false);
              break;
     }
  }
}

You can see the web site here: http://genepoulin.net/geo/index.html
Obviously the events are not triggered.

What am I missing here?

UPDATED and Resolved...

Corrected code (I was wrongly naming the switch cases):

function init(){
 var map, kml_1Layer, kml_2Layer, kml_3Layer;
...
kml_1Layer.events.register("visibilitychanged", kml_1Layer, layerChanged);
kml_2Layer.events.register("visibilitychanged", kml_2Layer, layerChanged);
kml_3Layer.events.register("visibilitychanged", kml_3Layer, layerChanged);
...
}

function layerChanged(layer) {
     if (layer.object.visibility) {
     switch (layer.object.name) {
        case "*given label of layer1*":
              kml_2Layer.setVisibility(false);
              kml_3Layer.setVisibility(false);
              break;
        case "*given label of layer2*":
              kml_1Layer.setVisibility(false);
              kml_3Layer.setVisibility(false);
              break;
        case "*given label of layer3*":
              kml_1Layer.setVisibility(false);
              kml_2Layer.setVisibility(false);
              break;
     }
  }
}

You can see the result on my page.


texts: How to change just one title when multiple titles are present in mxd?

texts:
I have an mxd with multiple titles and I'm trying to write a script that iterates a process that activates a layer, changes one of the titles, then exports it. I've got it all working except for changing the title. I can change the title but because there are multiple titles, it changes all of the titles to the same thing. Anyone know of a way to just target one of the titles within a mxd layout?

import arcpy

inputList = arcpy.GetParameterAsText(0)
bckGrndLyrs = arcpy.GetParameterAsText(1)
outputLocation = arcpy.GetParameterAsText(2)

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]

for lyrs in arcpy.mapping.ListLayers(mxd, "", df):
    if lyrs.name in inputList:
        lyrs.visible = False
    elif lyrs.name in bckGrndLyrs:
        lyrs.visible = True
    else:
        lyrs.visible = False

for inFC in arcpy.mapping.ListLayers(inputList, "", df):
    if inFC.name in inputList:
        arcpy.AddMessage("Activating " + inFC.name + " Layer")
        inFC.visible = True

        for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
          mxd.dataDrivenPages.currentPageID = pageNum
          for titleItem in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
            arcpy.AddMessage("...Changing title to " + inFC.name + " " + str(pageNum))
            titleItem.text = inFC.name
          arcpy.AddMessage("......Refreshing Dataframe")
          arcpy.RefreshTOC()
          arcpy.RefreshActiveView()
          arcpy.AddMessage(".........Exporting " + inFC.name + " Page " + str(pageNum) + " as PDF")
          arcpy.mapping.ExportToPDF(mxd, outputLocation + " " + inFC.name + " " + str(pageNum) + ".pdf")
        inFC.visible = False

del mxd


texts: Difficulty connecting QField to connect to PostGIS

texts:
I've been working on setting up a PostGIS server for QField for the past week, today I almost thought I had it working but now QField keeps throwing me one error no matter what I do:
"Connection to database failed could not translate host name 'My server IP address' to address: No address associated with host name"
I have the project working wonderfully in the QGIS environment. I'm not using the built in authentication, I'm using the pg_service.conf approach. I have all my layers set to 'no action'.
At one point I though I had the authentication setup right, then I tried getting Qfield to load a Postgres hosted project and stuff stopped working.
Does anyone have any ideas of a setting I might have screwed up or a step I'm missing?


texts: Limit the number of features rendered in a map

texts:
I have a layer of village points and I want to limit the number of features rendered to a maximum of 50 points. I tried using Layer Properties > Labels > Rendering > Limit number of features to be labeled to 50. But this obviously limits only the labels, and not the points themselves, like in the image below:

I also tried Layer Properties > Style > Control feature rendering order but it does not seem to do anything about the limiting of the points rendered in the map. How can I limit the number of features to be rendered, preferably using an expression/the Field calculator (so that I can dynamically set the same for other similar layers)?
I am using QGIS 2.18.


texts:
My requirement is to find a reliable source of Ocean water shapefiles, this is because OSM or Geofabrik does not provide this information.

I eventually landed on this page:
https://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-ocean/

By clicking on "Download Ocean" I was able to obtain the shapefile for the global Ocean. I've analysed it on QGIS, however, it turns out that, it is one big polygon.

My problems are:

1) I want to be able to refer to this in parts instead of one big polygon, I want to be able to cut out an area of interest of the ocean.
I am using QGIS, I have tried several clipping tutorials, but I must be doing something wrong, I keep getting the "0 features" error.
Please recommend, if possible with steps how I may achieve this with this specific shapefile. (Area of interest anywhere in the globe you may choose)

2) Is there another source for better segmented out Ocean water shapefile?

3) Is this a good source for water shapefile data?, if not, what do you recommend?


texts: Uninstall QGIS 2.18 from mac and reinstall. Remembers settings

texts:
Relates to Mac OS 10.14.

I have problems with QGIS 2.18. Plugins not working etc.
(latest problem is not being able to make a virtual raster catalogue of landsat 8 tiles. Error message :
"The process failed to start. Either the invoked program is missing, or you may have insufficient permissions to invoke the program.")

So I uninstalled QGIS by dragging QGIS.app from applications folder to trash and emptying trash.

Reinstalled QGIS 2.18.

Opened QGIS 2.18.

It remembered all the settings etc !

How can I delete QGIS 2.18 to start completely fresh ?
I imagine it's okay to keep the QGIS-OSX-2.18.28-1.dmg (this has all the packages required for QGIS to run)
I also have QGIS3 installed on the machine (also problems there) but I take it the two don't interfere, so I don't need to delete QGIS 3?

All help welcome.


texts: Coordinate transformation npm for X,Y and Z values

texts: Is there any direct node parameter libraries (npm) available to convert the X, Y, and Z values from one coordinate system to another? And are there any other methods for this conversion?


texts: Is it possible to get features from one layer within the bounding box of other layer?

texts:
I have 2 seperate layers say: water bodies layer, land layer from different sources but for the same location. as shown below:

No I want to use the bounding box of water bodies layer and get the featureID's of the 2 features in land layer.

I am using geoserver for testing. Is this possible using HTTP POST request? as my wfs doesnot support HTTP GET request. Please provide with any examples for better understanding.


texts: Field Calculator does not notice selection?

texts:
I have a problem with the Field Calculator in ArcGIS 10.

I selected specific points of a Shapefile I want to calculate, but FC is always calculating all fields. Unfortunately there is no checkbox in ArcGIS 10 anymore to tell FC to only calculate selected fields.

I discovered that it seems to be a problem with the Attribute Table since values can only be changed manually while being in the "Show selected records" mode. Otherwise my edits won't be even recognized (even actual typed in ones).

I tried both activated Editor and not activated Editor for FC. Just a few minutes ago everything worked out fine, now all Shapefiles can't be calculated according to the selection anymore. :/

Anyone an idea what went wrong?

I'm not aware of having changed any settings.

Restart of ArcMap also doesn't do it.


texts: Python script in Qgis for hundreds of texts replaces or gigant 'replace'+'case then else' combined expression?

texts:
I want to replace a text from one collumn into another text in a different collumn and base the condition on two different input collumn/fields.

Collumn 1: Integers (from 1 to 10)

Collumn 2: 100 different combinations of this text type - '1 and 10'

New collumn:

  When appears '10' on collumn1 and '1 and 10' on collumn2 replace with
  'YPC5' AND ( other example) when appears '1' substitute by 'YPC5' and when appears '10' substitute for 'YPC8', becoming: 'YPC5 and
  YPC8'. But, if the field1 changes for '9' then 'YPC5 and YPC7'.

Is there a python script for doing this in a easier way than use a gigant 'replace(' and 'case then else' combined expression?


texts:
I have converted a CAD file to GIS format using GIS (from .dgn to .shp)... In one layer I have polygons in which a smaller one is on top of the larger polygon. I can see both polygons however when I select the smaller polygon, the software automatically selects both the smaller polygon and the large one.
IS there anything I could do so that when I select the smaller one, the larger one doesnot get selected?

I don't want to select both and then unselect one of them.


texts:
I am loading data from a GPKG to a PostGIS database using GDAL and have noticed that it is very slow even with using the normal faster tricks. I am currently using this command
ogr2ogr -progress --config PG_USE_COPY YES -gt 65000 -f PostgreSQL PG:"" data.gpkg -nln output -nlt POLYGON -lco SPATIAL_INDEX=None

it has taken over 20 minutes to load in 65000 features.
Writing the other way is very quick so it seems to be the read from GPKG is bad.
Is there anything else I can try or do to improve the read rate?


texts: Custom metadata name in metadata.txt and plugins.xml

texts:
It is possible/recommended to add a new metadata name in metadata.txt and plugins.xml or is a closed format and we can have only the metadata name specified in Plugin metadata section ?

For example, I would like to add

<environment>SANDBOX</environment>

and the code will be:

<?xml version="1.0"?>
<plugins>
    <pyqgis_plugin name="Pluginname" version="plugin.4.0.400">
        <description>Plugin Description here</description>
        <homepage>http://www.plugi-home-page.uk</homepage>
        <qgis_minimum_version>2.14</qgis_minimum_version>
        <file_name>Plugin.zip</file_name>
        <author_name>Author name</author_name>
        <download_url>http://repository-url/QgisPlugins/plugin.zip</download_url>
        <environment>SANDBOX</environment>
    </pyqgis_plugin>
</plugins>

The startup.py will do an automatic upgrade of the plugin, if a new version is available or if the user changed the current repository (in current situation we are using multiple repositories, one for each environment: sanbox, prod... etc), so. a new tag to tell us for which environment is the available plugin, will be very useful.


texts: OpenLayers 3 - submit SLD text in WMS call

texts:
Earlier versions of OpenLayers allowed for the submission of the text of an SLD in a call to  the server side WMS.

This describes the procedure, and cautions that, if your SLD is too long, you might have to use POST rather than GET.

Documentation on how to do this with OpenLayers 3 seems to be hard to find.  Is it still possible to do this?

I did find that you can specify an SLD defined & stored within GeoServer by way of the STYLES property of the params argument to, for example, the ol.source.TileWMS.

I need to dynamically create the text of an SLD on the client side (within JavaScript) and submit this in my calls to get tiled maps from GeoServer.

Prior to OL 3, you could submit your SLD text with the SLD_BODY: parameter, but this does not seem to work, at least in the ol.source.TileWMS constructor.  I guess the crux of my problem is that the params object is not very well documented.


texts: One needs to be run from server and other runs itself, Why?

texts:
There are two types of project given in samples of developers.arcgis.com

https://developers.arcgis.com/javascript/latest/sample-code/source-code.zip
https://developers.arcgis.com/javascript/latest/sample-code/get-started-mapview/index.html

When I try to run from the latter(No. 2), it gives me GET Err:404_file_not_found but the former(No. 1) runs smoothly. But No. 2 works only when it is in the server (in my case I'm using XAMPP).

Why does one need to be put in the server to work and the other not?

I'm new to this field but have decent knowledge in javascript.


texts: How to display an arrow pointing to a specific bearing in QGIS

texts:
I have a number of point features which contain an attribute I call "bearing" and another one called "slope". Each point has a bearing between 0 to 359 degrees. And the slope is either "+15" or "-15", etc, such as +15% or -15% depending on the steepest slope of the ground.

What I am wanting to do is display an arrow on the map with the arrow pointing in the bearing direction. With a negative slope% pointing outward from the point and a positive slope% pointing inward to the point.

I'm using this in a forestry context with cruise plots. I'm fairly certain that ArcGIS can do this (As I've seen it on maps) but I'm hoping that QGIS can as well.


texts: store QgsVectorLayer.getFeatures() result in order to speed up subsequent feature requests

texts:
I'm doing a lot (really a lot) of feature requests on a vector layer. These requests are capsuled in iterations and functions, so no way of reducing the total number of requests.

Although I consider the #1 tip concerning QgsVectorLayer.getFeatures() - use a specific request whenever possible (Why use QgsFeatureRequest() when iterating?) - I end up with low performance and cProfile tells me the bottleneck is still the feature request.

Btw, using layer.getFeature(fid) instead of

request = QgsFeatureRequest().setFilterFid(fid)
feature = next(layer.getFeatures(request))

doesn't make any difference.

How do I speed up my feature requests, when I know each feature is called lots of times?


texts: Arcpy.mapping Select by Location problem

texts:
I need to adjust a premade .mxd many times, so I'd like to speed it up with this script. Everything is working except the Select by Location part.
The query reduce the Address layer to a single point and I want to use this point to select the layer. The mxd has all the county layers oppened but the input will change in the future.
I'm getting the error:ERROR 000368: Invalid input data. ERROR 000732: Selecting Features: Dataset Address does not exist or is not supported
If I run this script in ArcGIS with mxd set as "Current", then everything works.
I want to run the script without opening the actual mxd.

import os
import shutil
import glob, os

SID = "12234"
SIDtext = "df-" + str(SID)
address1 = "1233 Antler Dr"
address2 = "Summerton"
address3 = "23332"
address = address1 + "\n" +  address2 + "\nSC " + address3
county = "Clarendon"
owner = "Smith"

test = "C:\\Users\\test\\test.mxd"

mxd = arcpy.mapping.MapDocument(test)
for df in arcpy.mapping.ListDataFrames(mxd):
    #adjusting 1st dataframe
    if df.name == "Overview Map":
        for lyr in arcpy.mapping.ListLayers(mxd):
            if lyr.name == "Address":
                lyr.definitionQuery = "Case_ID = '" + str(SID) + "'"
    #adjusting 2nd dataframe
    else:

        df2 = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
        for lyr in arcpy.mapping.ListLayers(mxd, "Address", df2):
            if lyr.name == "Address":
                lyr.definitionQuery = "Case_ID = '" + str(SID) + "'"
            df2.panToExtent(lyr.getSelectedExtent())

        arcpy.SelectLayerByLocation_management(county, "INTERSECT", "Address")
        arcpy.MakeFeatureLayer_management(county, "Property Parcel")
        for lyr in arcpy.mapping.ListLayers(mxd):
            if lyr.name == "Property Parcel":
                lyr.visible=True
                symbologyLayer = "Terra selection"
                arcpy.ApplySymbologyFromLayer_management (lyr, symbologyLayer)

        arcpy.SelectLayerByAttribute_management(county, "CLEAR_SELECTION")
        arcpy.RefreshActiveView()

mxd.save()
del mxd

Error message:

Runtime error  Traceback (most recent call last):   File "", line 33, in    File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 7320, in SelectLayerByLocation     raise e ExecuteError: ERROR 000368: Invalid input data. ERROR 000732: Selecting Features: Dataset Address does not exist or is not supported


texts: Making paths and roads work together with ArcGIS Network Analyst?

texts:
Working in Arc GIS 10 .0 service pack 5, network analyst extension, Productivity Suite 2.1 service pack 7, on Windows 7.

I am trying to get the Ordnance survey ITN urban paths layer to work with OS ITN roads to model pedestrian journeys. I've had a ton of problems doing this mainly related to the sofware that unpacks the OS raw data, but I think (?) it finally unpacked it OK and my problem is in Network Analyst in Arc.

I've got threads on a bunch of forums and I think I'm inching towards a soloution. I need nearest facility and service areas that are based on connections across both paths and roads in a network. I had a network that would do solves but only on paths or roads but no connectivity between them (the connectivity properties were fine though).

From other advice I got I found there were two problems here, firstly a bug in Productivty Suite which the service pack (and a reinstall) got around.

Secondly I'd not integrated the two layers - so I did this using the integrate tool without crashing or errors, although I'm not sure I've done it right - I just integrated the path links and road links - not sure if I need to select all of the nodes as well and integrate the lot?

Assuming I have done this correctly, I now try a simple route solve to test the data that have been unpacked and integrated and I get an error message that says it is "unable to create a network anlalysis layer because the network data set has no cost attributes".  Not sure where my cost attributes went or what to do about this. Incidentally I'm not doing anything very sophisticated with costs, I simply want a bunch of nearest facilities and service areas based on walk distances in metres - no parameters use of hierarchies or anything.

I've got the technical paper from the OS on using Connecting Links in Urban Paths.

Has anyone out there made these data work together?


texts: NNJoin in QGIS does not return any distance

texts:
I have two layers, one with ~100,000 points, and another with ~20 polygons.
I'm trying to find:
a. which polygon a point belongs to
b. the distance from the point to that polygon
I understand that NNJoin should be able to do both of these. However, it only gives me the answer for (a). The "distance" column - which should contain the answer to (b) - is totally filled with 0s.
When I run NNJoin, I get a warning message saying "Geographic CRS used for the join layer - distances will be in decimal degrees" - but I get 0s, not decimal degrees. I've checked that I'm using the same CRS throughout my project (EPSG:4326), so I'm stumped (I get perfect distance values, in metres, using the "distance to nearest hub" tool). Does anyone know how to fix this?


texts: Is there a way to turn on labels with ModelBuilder?

texts:
I have a model that creates a feature class and adds it to the map. The properties of the feature class reference a layer file for the symbology. However, I can't seem to find a good way to pass the labels as well so that the feature class is added with the labels visible.

Is there a way to do this?


texts: Reading Esri file geodatabase (*.gdb) using DotSpatial?

texts:
I am a big user of DotSpatial and have even contributed some source code. However, I am having some roadblocks in trying to figure out if DotSpatial can read a file geodatabase. I have posted on the official forum but don't think its very active right now.

Has anyone been able to do this with DotSpatial?


texts: BoundingBox of Layer seems to be unnaturally large

texts:
I retrieve the bbox of my main Layer via GetCapabilities.
Now, if I instantiate my map like this, I would expect to have the map centered within the bbox.
map = new OpenLayers.Map({
    div : "map",
    projection :  new OpenLayers.Projection("EPSG:31466");
    units : "m",
    maxResolution : 2,
    maxExtent : new OpenLayers.Bounds(bbox)
});

However, when I do it like this, the map is displayed in the upper right corner.
Even stranger, when I put markers on the edges of the bounding box, the upper right corner is very far away and the lower left is not visible at all. (See picture, the red markers are the markers on the edge of a hard-coded bounding box that I use and that I want to replace by a Server-generated Version, the green marker is the upper right of my layers bounding box according to GetCapabilities).

EDIT:
The bounding box that getCapabilties sends:
[360000, 5770215.8928236235, 365000, 5770804.068148668]
The bounding box that is hard-coded (red markers): [363060,5770362,363496,5771118]


texts: Spatial search to find featureclasses in SDE

texts:
Is there a way to find all featureclasses within SDE that fall within a spatial search e.g a user draws a rectangle on screen over a map of the world say.  (I would like to do this progrmammatically, the native search tool in ArcCatalog is horrible).

I'm thinking you may be able to do something with ArcSDE system tables to find things out or maybe go down to the level of creating some custom code in C/Java SDE API.

So this is like doing a spatial search WITHIN a featureclass to find features - but a step further up in thinking.

Many thanks.


texts:
I've Mapinfo .tab files (including .dat .id .map) I'd like to read the attributes of the files then may store in a DB or write to a csv (not yet decided).

As far as I can tell, using  MapInfo Pro will give me those attributes.

I'm just wondering is there any open source that will do the same?

I'm new to GIS.


texts: How do you count the number of pixels that have overlapping data for a single land area in Google Earth Engine? Count No-null values?

texts:
I have an Image Collection of Images that have pixels with values of either 1000, 1000, or 1. I want to know which pixels, over the years have a value of '1000' 75% of the time, which pixels have a value of '100', 75% of the time, and which pixels have a value of '1' 75% of the time. Time being years (each image in the collection represents a single year). However, some Images do not have pixels for that year (so we have null values I assume). So, is it possible to find the 75% without knowing how many years of actual data (does GEE just ignore the null values if I were to run a calculation on the Image Collection), or... do I need to answer the following question?
How do you find out for every pixel (land area), how many years of data an Image Collection stacked together has? Some years do not have data (i.e., a pixel may have data 4/5 years in the Image Collection, another pixel may have data 3/5 years in the Image Collection). How do you do this in Google Earth Engine and then store it as a variable I can do calculations on? Or how do you count the non-null pixels and return them as an image for calculations?


texts: Best Method to Replicate ArcPy Toolbar in PyQGIS

texts:
I am trying to build the equivalent of an ArcMap Toolbar in PyQGIS. My toolbar has interactive functionality. For example, in the toolbar below, I cannot access any edit functions until I have selected a certain raster id. Each button has a different purpose and I still need to be able to draw polygons on the map canvas. I need a similar functionality in my QGIS version.

I also need to be able to override some standard QGIS functionality such as some forms when saving layer edits ( I understand that I can do this using slots and signals)

Currently I have come across 2 different methods of trying to replicate this in QGIS:

Create a plugin for each button using the QGIS Plugin Builder and Qt Designer and then group my plugins into a toolbar using a Custom Toolbar Plugin. Although this is what I need aesthetically, it doesn't allow interactivity between my plugins(???)
Use PyQt to build a toolbar from scratch. Although I think this will give me the functionality I need, it seems to be far more convoluted and as soon as I click off the toolbar it disappears- I need it to be present at all times.

Can anyone advise on what the best method to approach this problem is? I.e. what is the best method to get full interactivity between the different buttons/plugins and yet group them together as a toolbar?

I am using QGIS 3.4 and PyQt5


texts: QGIS not responding when committing changes to a layer

texts:
I run QGIS version 3.10.6 on my laptop, which runs Windows 10. The program is a bit slow on startup, but the real issue comes when I calculate a field on a layer to add to my attribute table. When I exit the edit mode and save these changes, the program stops responding.
I have already tried to disable plugins, which hasn't changed the situation. This is a fresh install on Windows.


texts:
I need to develop a website to find the fastest road between two points (just like google map)
What kind of technologies I need to use? What are the components of this project (database, website, etc)? What else do I need to consider?


texts: Making contour lines out of points with different heights using QGIS?

texts:
I have different points in a polygon. The points have the attribute "depth" and either have 0m, 2m or 4m depth. My task is now to create some sort of contour lines with the points with the same depth so that it is possible to see how deep the polygon is.

How can i do this?

I work with QGIS 2.4.


texts: Difference between Vincenty and great-circle distance calculations?

texts:
Python's geopy package features two distance measurements techniques: Great Circle and Vincenty's formulae.

>>> from geopy.distance import great_circle
>>> from geopy.distance import vincenty
>>> p1 = (31.8300167,35.0662833) # (lat, lon) - https://goo.gl/maps/TQwDd
>>> p2 = (31.8300000,35.0708167) # (lat, lon) - https://goo.gl/maps/lHrrg
>>> vincenty(p1, p2).meters
429.16765838976664
>>> great_circle(p3, p4).meters
428.4088367903001

What is the difference? Which distance measurement is preferred?


texts: Drag and drop features of a WMS service using OpenLayers

texts: I have a WMS service, which is a point layer. Is it possible to drag and drop features of a this service layer through OpenLayers?


texts: ECEF coordinates to geodetic coordinates: Where is error in calculation?

texts:
I know, there are several questions to this topic, but I am just searching for the error in my own equations. I have no deeper geodetic background, but I am trying to understand these formulas instead just copy them.

So, I have a XYZ coordinate and want to get a geodetic coordinate. This is what I tried:

tan(longitude) = Y/X

For the latitude I tried this approach. On a sphere you could calculate it that way:

diameter_xy = sqrt(y^2 + x^2); // => d
tan(latitude) = z / d

Now, I have no sphere but an ellipsoid. So the latitude I already calculated is the same as the geocentric latitude of an ellipsoid, right? (let's call it p)

I tried to transform this into the geodetic latitude (q). The book "Geometric Geodesy I" from R. Rapp (p. 25, formula 3.62) creates this relation between geocentric and geodetic latitudes:

tan(p) = (1 - e^2) * tan(q)

So:

tan(q) = tan(p) / (1 - e^2)
       = (z / d) / (1 - e^2)                // taking the spherical formula above
       = (z / sqrt(y^2 + x^2)) / (1 - e^2)
       = z / ((1 - e^2) * sqrt(y^2 + x^2))

Well, now I don't come to an expected result. When I have a look a several solutions, the formula differs:

For example: https://gis.stackexchange.com/a/292635 The initial equation for the latitude seems to be the same as mine (beside the point, that atan2 is used instead of atan, which is also not that clear to me). But later the height h becomes part of the equation. It absolutely not clear to my, why.

Well, I don't understand were is error and where do the real solutions (e.g. Ferrari or so) get that complicated.

Edit:

In my case the values result in a greater height. The XYZ point is not on the surface. But since I am not interested in the height above ground, but only in the angles, I thought, this wouldn't make any difference. Maybe I am wrong? Does my calculation only work for h == 0m above ground? If so, how do I recognize the height and how do I eliminate this?

My use case is:

X = 7000000m
Y = 2000000m
Z = 10000000m
e = 0.00669438 // WGS84

So, I am getting:

lat: 53.95°
lon: 15.95°

The lat real value is something like:

lat: 54.04° (h: 6,005 km)


texts: How to debug GeoServer HTTP response code: 400?

texts:
We are using opengeo 2.4.3, and trying to open some layers with GeoServer through WMS. The different layers in the WMS file are detected, but when doing 'layer preview', we get the following error:

'org.geoserver.wms.map.RenderedImageMapOutputFormat.directRasterRender(RenderedImageMapOutputFormat.java:871)
at org.geoserver.wms.map.RenderedImageMapOutputFormat.produceMap(RenderedImageMapOutputFormat.java:290)
... 88 more
Caused by: java.io.IOException: Server returned HTTP response code: 400 for URL:
http://localhost:8080/thredds/wms/conakry.nc?
SERVICE=WMS&
LAYERS=tz&
FORMAT=image%252Fpng&
TRANSPARENT=FALSE&
HEIGHT=512&
BGCOLOR=0xFFFFFF&
REQUEST=GetMap&
WIDTH=512&
BBOX=-21.017490386963,3.9957752227786,-9.0176734924312,15.99559211731&
STYLES=&
SRS=CRS:84&
VERSION=1.3.0'

The thing is that we have previously done this with the GeoServer 2.1.1 version and everything worked fine. We have checked the URL causing the problem, and pasting the URL into a web-browser and changing the 'SRS' bit to CRS bit seems to be working somehow...

Can someone please help us on this issue?

Otherwise, does somebody know how to use another version of GeoServer in the OpenGeo Suite


texts: How to set the properties of L.Icon.Default.prototype.options? I want to change the size of L.Icon.Default.


texts: Getting list of selected features in ArcGIS for Desktop using Python code?

texts:
I have a question regarding selections in ArcGIS for Desktop. Assumed I have one layer in ArcMap and I have selected two of five features.

Is it possible to get a list of all selected features by using Python?

It would be fine if there is a way to get one special (or all) attribute(s) of the selected features stored in a list that can be written into a txt file.

Is it possible to do this in ArcGIS for Desktop?


texts: PostGIS - create polygons from clusters of a fix number of points

texts:
I have a set of about 120,000 points, not regularly spread.

I would like to cluster them in order to create groups of about 3,000 points, regardless the distance. The final goal is to create polygons (their size is not important, I just want them to contain 3,000 points) out of bbox from these groups. The result would be like a grid with different sizes of cells, each cell containing 3,000 points. Do you know how could I achieve this ? I saw st_clusterKmeans and st_clusterDBSCAN functions, but I don't know how to fix the number of points to aggregate (there's a minpoints with st_clusterDBSCAN, but no "maxpoints").


texts: I can't find geotag and import photos plugin (QGIS plugin)

texts: I am doing a project on mapping tourist attractions and would like to attach photos and videos to point locations. I have read about the Geotag and Import Photos plugin but I can't find it in my plugin lists. I am using QGIS 2.8.2 Wien and I have downloaded the exiftool already.


texts: Why WMS works and WMTS throw out of bounds

texts:
I'm having some troubles getting wmts service to work. This is what I have done so far:

I have imported OpenStreet Map info for Spain using osm2pgsql into a postgis database. I get Spain info from Geofabrik
I have published planet_osm_line with this bounds:
I have enabled tile caching in the layer, with this config:

Then from layer preview menu, I can see values with no problem using open layer. You can check here

However, from tiles cached menu, if I try to open EPSG 900913 view, only errors are thrown in the chrome console, with an out of bounds message. You can check it here

    <ExceptionReport version="1.1.0" xmlns="http://www.opengis.net/ows/1.1"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://geowebcache.org/schema/ows/1.1.0/owsExceptionReport.xsd">
  <Exception exceptionCode="TileOutOfRange" locator="TILEROW">
    <ExceptionText>Row 131073 is out of range, min: 131071 max:131072</ExceptionText>
  </Exception>
</ExceptionReport>

Any ideas of whats going on? Why WMS service works and WMTS do not? What have I  to configure to get it work?


texts: Represent Ternary Legend for IMG in ArcGIS Desktop?

texts:
I've got a composite image built up from Red Green and Blue composite colours. Each channel represents the differing concentration of a different variable in each location - with 0 being the lowest concentration and 255 the highest.
I'm bringing each of these unique values together on the composite image to form a colour based on the RGB values.

That bit is straight forward - although I'm trying to find a useful way in which to represent the colour scheme derived as a Legend. I've tried various RGB triangles and cubes, although nothing seems to fit quite right.
I was wondering if anyone had attempted this before / could point me in the direction of some Creative Commons images that could be used to this end?
Heres a screen shot of what I'm trying to represent:


texts: How to write a string based rule, e.g. If "Name" contains 'NRG'

texts:
I have created a field called 'Name' in QGIS and I wish to create a rule-based style for a gridded vector. The grids are either 'NEP1' or 'NEP1NRG' etc. with values to 5000+. I want a rule to colour the grid polygons with 'NRG' in the 'Name' field one way, and those without this in the name another.
E.g. If "Name" contains 'NRG' then adhere to style rule, if not then another.

I have tried variations of LIKE but not sure this is appropriate.


texts: Select feature one by one and save it as image

texts:
I have used following code for selecting the features and save it as image.

layer = self.iface.activeLayer();
layer.setSelectedFeatures([])
count =layer.featureCount()
lstChalta = []
lstFeatures = processing.features(layer)
for elem in lstFeatures:
   geom = elem.geometry()
   attr = elem.attributes()
   layer.setSelectedFeatures([elem.id()])
   self.iface.mapCanvas().zoomToSelected( layer )
   box = layer.boundingBoxOfSelected()
   self.iface.mapCanvas().setExtent(box)
   self.iface.mapCanvas().refresh()
   self.iface.mapCanvas().refresh()
   cno= str(attr[layer.fieldNameIndex('cno')])
   path = 'D:\Test'+ str(cno) + '.png'
   self.iface.mapCanvas().saveAsImage(path)
   #QMessageBox.information(None, "DEBUG:", str(test))

Whenever I run the script, it is selecting only the last feature. However, if I uncomment message box line it is working fine and selecting each and every feature one by one. In addition to that it is saving the image properly.

I want to work the script without showing message.


texts: Aggregate attribute based on symbology classes in QGIS

texts:
I was wondering if there is an existing tool/way to aggregate an attribute based on the active symbology of a layer.
For instance, if I have a point layer with population and distance from the city center for each point and I create a specific symbology for the distance (let's say 0-1km, 1-5km 5++km) to get the total number of people within each class.
I would like to avoid creating dummy reclassified attributes if possible.


texts: Running Python script in Task Scheduler-- Script will not run

texts:
I have a script which I wish to run on a regular monthly basis without me needing to open and run the it.  I have referred to every web page I can find on Windows Task Scheduler and followed their instructions.  However, running the script from Task Scheduler will only open the script in Python and not run it.

In the Task Scheduler Run Window, I have tried entering both the path to the script alone and the path to the script with the path to my Python program.  I've also tried entering the path to my Python program as the script with my script as an argument as suggested in this blog post. My script is hard-coded, so I don't think I need to add any arguments.

Am I missing something obvious?

I was experiencing issues because I was using the path to the Python program and not the Python command line.  I.e., I was using "C:\Python27\ArcGIS10.1\Lib\site-packages\pythonwin\Pythonwin.exe" when I needed to be using "C:\Python27\ArcGIS10.1\python.exe".


texts: Compare two tables. Show values only in one table

texts:
I have two tables. Each table has a column with a reference number. Most of the records in these tables are the same but there are some which are in one table but not in the other.

How can I compare these columns and pick out the records which are only in one table?

I could write a MapBasic program that loops through each record and looks it up in the other table but I'm wondering if there is a simpler SQL-like way of doing this.


texts: Unable to edit geometry into database using PostGIS connection in Civil 3D

texts: I have prepared a set of tables (in PostgreSQL) in database with unique geometry types (point, line, polygons) for survey work. I am using Autodesk Civil 3D to edit the geometries in the tables but they are not getting saved in the database (which is feature source). However, I am able to create a geometry and edit it once in Civil 3D and save in the database (using check in and check out) but i am unable to make any further changes. I am working by connecting to the PostgreSQL connections in planning and analysis of Civil 3D workspace. I was wondering if anybody has faced the same issues in a similar type of workflow.


texts:
I've been using QGIS 2.6 and whenever I try to install some plugin through plugin manager, it returns an error "Plugin installation failed:SSL handshake failed".
I have always had to manually install plugins since then.

Can anyone suggest me solution to this ?


texts:
I'm new to ArcScene so maybe I'm doing something wrong, but I can't see any colors in my ArcScene views. I've attached 3 screenshots, one from ArcMap showing DEM and hillshade working just fine, and 2 from ArcScene that show a lot of grayscale rendering in both raster (DEM) and vector layers (streets, buildings) despite changing the symbology in preferences.

Any ideas?


texts: Warp vector WGS84 polygons to an orthorectified sub image

texts:
I'm using GDAL Python bindings (and gdal command line tools via subprocessing) and I'm having trouble warping WGS84 polygons to an orthorectified subimage. I'm fairly new to GIS, so please excuse any incorrect terminology. I'll try to explain in as much detail as possible.
First, here is what I have:

A large geotiff with RPC metadata.
A WGS84 region of interest.
Multiple WGS84 polygons inside that region of interest.

What I have been able to do

Use gdalwarp to extract an the orthorectified ROI from the large geotiff. I did this via gdalwap:

gdalwarp
-te {xmin} {ymin} {xmax} {ymax}
-te_srs epsg:4326
-s_srs epsg:4326
-rpc -et 0
-to RPC_DEM={dem_fpath}
-overwrite
{orig_geotiff} {roi_geotiff}

Where xmin/ymin/xmax/ymax are the lat/lon bounds of the ROI, dem_fpath is the path to the corresponding DEM tif, orig_geotiff is the large source geotiff and roi_geotiff is where the new much smaller orthorectified image is written.

I can warp the WGS84 polygons into pixelspace of the original orig_geotiff large geotiff.

import gdal
src_ref = gdal.Open(orig_geotiff)
tr = gdal.Transformer(src_ref, None, ['METHOD=RPC', 'RPC_PIXEL_ERROR_THRESHOLD=0.00' f'RPC_DEM={dem_fpath}'])
orig_pixel_polys = []
for poly in geo_poly_list:
    pts_in = poly.data['exterior'].data
    pts_out = []
    lats, lons = pts_in.T
    # I can use the DEM to find a better
    # initialization for this but I'm keeping it
    # simple in this example
    elevs = [0] * len(pts_in)
    for x_in, y_in, z_in in zip(lons, lats, elevs):
        bDstToSrc = 1
        success, pnt = tr.TransformPoint(bDstToSrc, x_in, y_in, z_in)
        x_out, y_out, _ = pnt
        pts_out.append((x_out, y_out))
    pts_out = np.array(pts_out)
    orig_pixel_polys.append(pts_out)

What I need help with
But what I'm having an incredibly hard time doing is trying to either (1) warp those pixel-polygon points in the original image to the orthorectified image or (2) warp the WGS84-polygon points directly into the orthorectified image.
In the first case, I'm not sure how to find whatever the transform gdalwarp is using and apply it to a vector dataset.
In the second case, the file written by gdalwarp does not recompute new RPCs for the orthorectified subimage; it only brings over the affine geotransform. Unfortunately, this affine transform does not do a good job of mapping the world coordinate system to the pixel coordinates. If I look at the alignment between the data and the pixel coordinates in the original image it looks fine (because it used the RPCs to do map WGS84 to original pixel coordinates) but when I try warping using the affine geotransform in the orthorectified image (or in the original image for that matter) they do not align well.
Summary
In summary I'm looking for an answer for one of the following:

How do I take pixel coordinates in an original image and map them to the corresponding pixel coordinates in a a cropped subimage that was orthorectified with RPC metadata?

How do I take WGS84 coordinates and map them into a pixel coordinates in a cropped subimage that was orthorectified with RPC metadata?

A method for either of these would answer my one question of: "How do I warp vector WGS84 polygons to an orthorectified sub image".


texts: Calculating a baseline shoreline from multiple surveys

texts:
We are trying to calculate a average shoreline position from a collection of shorelines dating from the 1880s to the present, the purple lines shown on the map below:

We want to create some code that averages them to create a baseline shoreline, something like the hand-drawn black line on the map.  Ideally the resulting line should fall in the middle of the shorelines and also cross inlets.

We were thinking to just take offshore transects every so often, compute the mean, and then spline fit the resulting points, but it seems there are lots of things to consider here that could cause problems.

Does this approach sound reasonable?
Or is there a better way to do this, ideally in python?


texts: GDAL Calc using different raster grid sizes

texts:
I have several GeoTIFFs (.tif) that I am conducting raster calculations on using gdal_calc.py, all three are of different sizes. These rasters can be best illustrated with the following image:

All three rasters domains overlap across the northern extent, however the blue raster is of full size and covers the entire area of interest.
I'd like to perform a simple raster calculation (say gdal_calc.py -A redgrid.tif -B greengrid.tif -C bluegrid.tif --calc="(A+B+C)/3"), and generate an output raster in which I will get resultant values across the extent of the largest raster available, or the blue domain. Currently, this calculation will only yield an output raster with the smallest extent, or the red domain.
I've considered filling the other two grids with np.nan to match the extent of the blue domain, though I believe that this division would return nan in areas where the three grids do not overlap.


texts:
I have a map that is using OpenLayers and GeoServer. I have a base layer and a WMS layer with many features. I was using OpenLayers.Control.Click to get the coordinates from the map and show a popup about the closest feature on the WMS layer.

This worked some of the time but it didn't always bring up the popup for the selected point. When it did work, though, it worked on desktop and mobile. So I changed to OpenLayers.Control.GetFeature. This works great on the desktop but doesn't respond at all on mobile.

Can I get it to work on mobile?


texts: How to embed a mapview map in a Jupyter Notebook with IRkernel?

texts:
I'm looking to integrate some awesome mapview interactive maps in a Jupyter Notebook with the IRkernel.

I would like to know if:

Is it possible? How can I achieve that?
Is there any workaround using Jupyter Notebook with IPython kernel using cell magics %%?

In a Jupyter Notebook cell, I tried:

# Load shapefile
shapefile <- sf::read_sf('../path/to/shapefile.shp')

# Plot mapview map
mapview::mapview(shapefile)

No luck, and no error message showed!


texts: How do I statistically compare two maps?

texts:
I need to statistically compare two maps in order to determine if the spatial distribution of the data is correlated or not.

I have two files in ArcGIS

XY location of foreign-owned firms
XY location of domestically-owned firms

I want to test if there is a correlation between the location of foreign and domestic firms.

I am comparing points within the same map and within the same period (I have 7 time periods, a dynamic analysis would be nice, but too complicated for now)

How might I go about this in ArcGIS?

From my own research, I can see that there is no universal way to do this.

Some suggestions I have seen online are Mantel tests and Cross-covariance analysis

Or, Using the Geoda or Crimestat softwares (I have tried both, but I am still at a dead-end).

I just want to calculate what the correlation between the two maps is, I know there is many ways to do it, but I can't find the correct tools to do even just one of them.


texts: Points in Polygon Count: Error with arcpy.selectLayerByLocation_management

texts:
I'm working on generating some multivariate hexbins and am running into a problem with the initial step of counting the number of points (for the first variable) that are within each hex polygon. Here is the code that I have generated.

import arcpy

arcpy.env.overwriteOutput = True
arcpy.env.workspace =  "d:/projects/scf/shapefiles"

hexBins = "d:/projects/scf/shapefiles/hexbins.shp"
points = "d:/projects/scf/shapefiles/points.shp"

bid = "BIN_ID"
updateField1 = "UCR"

with arcpy.da.UpdateCursor(hexBins, (bid, updateField1)) as binRows:

        arcpy.MakeFeatureLayer_management(points, "pointLayer")

        for bin in binRows:

            binIDstring = bin[0]
            queryString = '"' + bid + '" = "' + str(binIDstring) + '"'

            arcpy.MakeFeatureLayer_management(hexBins,"CurrentBinLayer",queryString)

            selPts = arcpy.SelectLayerByLocation_management("pointLayer",'WITHIN',"CurrentBinLayer")

            getPointCount = arcpy.GetCount_management(selPts)
            binPointCount = int(getPoimtCount.getOutput(0))

            bin[1] = binPointCount
            binRows.updateRow(bin)

            arcpy.Delete_management("CurrentBinLayer")

        arcpy.Delete_management("pointLayer")

When I run this code, I get an error at the SelectLayerByLocation line.

  ExecuteError: ERROR 999999: Error executing function.
  A column was specified that does not exist.
  A column was specified that does not exist.
  Failed to execute (SelectLayerByLocation).

However, after running this code, I've double checked that the MakeFeatureLayers were valid (they are) and then I've successfully run the line

selPts = arcpy.SelectLayerByLocation_management("pointLayer",'WITHIN',"CurrentBinLayer")

I've tried moving the folder to different locations (testing for file permissions issues) but to no avail.

Amy thoughts?

Thanks.


texts: Loop through table, look for word in text field, set boolean true if word exists

texts:
I want to create a function that loops through a table and checks if a specific word exists in a text field. If it does, the function should set a Boolean field true. This is what I've come up with. I guess my execute command is wrong?

CREATE OR REPLACE FUNCTION insertCarBooleans()
Returns void
AS $$
DECLARE
    location RECORD;
BEGIN
    FOR location IN SELECT * FROM Ownlocations ORDER BY id LOOP
    EXECUTE
        'UPDATE Ownlocations SET Ford = true
        WHERE
        cars LIKE 'Ford';
    END LOOP;
RETURN;
END;
$$
LANGUAGE plpgsql;

SELECT insertCarBooleans()


texts: How to increase the number of features returned by an ArcGIS Online feature service

texts:
When uploading a feature layer to ArcGIS Online, and then accessing it via ArcGIS Online's REST API, there is a limit of 1000 features returned.

Is it possible to increase the limit and return more than 1000 features?


texts: Dynamically editing items in combobox of Python AddIn?

texts:
I'm trying to dynamically edit the list of items in a combobox I made using python add in for a toolbar I'm working on. The list is calculated by another button. Here is the relevant code:

class SetLayer(object):
    """Implementation for stuff.button2 (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        self.a = pythonaddins.GetSelectedCatalogWindowPath()
        self.a = self.a.encode('mbcs')#pretty sure this helps with string encoding
        print self.a
        self.b = os.path.split((self.a))
        self.c = self.b[0]
        self.d = os.path.split(self.c)
        self.e = (self.b[1]) #feature class
        self.f = (self.d[1]) #feature dataset
        self.g = (self.d[0]) #workspace
        print "Feature class: %r" % self.e
        print "Feature dataset: %r" % self.f
        print "Workspace: %r" % self.g
        connProp = arcpy.Describe(self.g).connectionProperties
        self.h = connProp.instance.split(":", 2)
        self.k = self.h[2].encode('mbcs')#instanceX
        self.i = connProp.database.encode('mbcs')#databaseX
        print "Instance is %r" % self.k
        print "Database is %r" % self.i
        print "Availible Versions:"
        self.List = arcpy.da.ListVersions(self.g)
        self.nameList = []
        #This creates the list that I want the combobox to use
        for version in self.List:
            self.nameList.append(version.name)

        print self.nameList
        print "Tool finished"
        pass

#working on editing this to provide dropdown
class setVersion(object):
    """Implementation for stuff.combobox2 (ComboBox)"""
    def __init__(self):
        #self.items = []
        self.editable = True
        self.enabled = True
        self.width = 'WWWWWW'*3
        self.dropdownWidth = 'WWWWW'*3
    def onSelChange(self, selection):
        self.items.append(button2.nameList)
        self.version = selection#change this
        print self.version
        pass

How do I get the items to dynamically update? What did I do wrong?


texts: Overlaying user raster datasets on nasa world wind

texts:
Currently I am developing a desktop application with Nasa world wind JDK for 3d visualization of my data sets, but when i overlay my raster data sets those are vertically exaggerating with the default elevation layer.Can anyone help me how i can visualize my DEM data(SRTM) on nasa world wind in 3d with its own vertical exaggeration.
Thanks in advance


texts: I've established how to run a Least Cost Path analysis in ArcGIS, but I'm wondering if you can run an analysis to find the hardest path between two points?


texts: R Reprojected raster stack loses band information

texts:
I'm attempting to reproject a multiband raster image in R to Orthographic for export
https://www.naturalearthdata.com/downloads/50m-raster-data/50m-cross-blend-hypso/

However when I project the raster I get the following error, and the result although reprojected is not the bands that I started out with

Error in if (value[1] != nrow(x) | value[2] != ncol(x)) { :
  missing value where TRUE/FALSE needed
In addition: Warning message:
In `dim<-`(`*tmp*`, value = c(nr, nc)) :
  NAs introduced by coercion to integer range

This is how I load the raster

library(raster)

# Create RasterStack object
r <- stack('HYP_50M_SR/HYP_50M_SR.tif')

In the console this returns

coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

        HYP_50M_SR.1 HYP_50M_SR.2 HYP_50M_SR.3
Min.              64           94           88
1st Qu.          199          199          193
Median           250          250          250
3rd Qu.          251          251          251
Max.             255          255          255
NA's               0            0            0

Then I reproject

# Define the Proj.4 spatial reference
sr <- "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"

# Project Raster
projected_raster <- projectRaster(r, crs = sr)

In the console this returns

coord. ref. : +proj=ortho +lat_0=-10 +lon_0=-60 +ellps=WGS84

      HYP_50M_SR
Min.    7.105387e+01
1st Qu. 2.027389e+02
Median  2.502865e+02
3rd Qu. 2.507040e+02
Max.    2.550000e+02
NA's    4.880675e+06

How can I reproject a multiband raster while keeping the band information?

Ive been using this for reference https://datacarpentry.org/r-raster-vector-geospatial/05-raster-multi-band-in-r/


texts: Importing raster dataset in ArcSDE/PostgreSQL/Postgis database useable in QGIS and ArcGIS?

texts:
I want to import a raster dataset in my database (PostgreSQL 9.4, Postgis 2.1), the catch is that I also want the raster to be viewed in QGIS via a Postgis connection.

Importing feature classes in my database, I can simply use PG_GEOMETRY to make these readable in QGIS, but how does it work with a raster dataset?

Simply importing a raster via the ArcCatalog does not work.


texts: 'pyrcc4' is not recognized as an internal or external command?

texts:
C:\Users\Department\.qgis2\python\plugins\np>pyrcc4 -py2 resources.qrc -o resources.py

  'pyrcc4' is not recognized as an internal or external command,

  operable program or batch file.

How do I solve this problem?


texts: Absolute scale for Geoserver's heatmap rendering transformation

texts:
By default the SLD to generate a heatmap in Geoserver takes numbers between 0 - 1 to define the color ramp. So this would be a relative scale, where 1 is assigned to the area with the highest value in relation to the rest of the dataset.

I have multiple layers and I want to visualize them on the same scale. Each layer is based on different data, I would like to visualize them on an absolute scale.

For example: I'm using the following color-ramp: green-yellow-red. I would like all area's with a value of 500 + to be red, with a value of 250 to be yellow and with a value of 1 to be green. This means I might have some layers without red, or yellow area's. I want this so I can visually compare the different layers to eachother.

Is there a way to do this, without writing my own rendering transformation in Java?


texts: Google Earth Engine - filter landsat by cloud cover over multiple polygons

texts:
I want to make a Landsat image collection that is totally cloud free over a SET of polygons. I have about 300 polygons (all within a single Landsat tile). I found this question, which is doing exactly what I want over a single polygon, but I can't figure out how to get it to work over a set of polygons (which are loaded as a feature class).

Here is what I tried to do to modify the solution from the question posted above, where imgcollection is my Landsat TOA Image Collection and fc is my feature class of polygons. I was trying to use map and clip together to select all the polygons in the feature class for each image (based on the answers in this related question).

  var combine = imgcollection.map(function(imgcollection) { return
    imgcollection.clip(fc); });

After this line, I tried to continue with the solution from the first question, but it doesn't work.
The error I get is as follows, which confuses me. I'm trying to use the variable "combine" as the geometry - am I misunderstanding this and actually using something else?

  ImageCollection (Error)
  Image.reduceRegion, argument 'geometry': Invalid type. Expected:
  Geometry. Actual: ImageCollection.

EDIT: full code

  var fc = ee.FeatureCollection('ft:omit');

  var imgcollection = ee.ImageCollection("LANDSAT/LC08/C01/T1_RT_TOA");

  var combine = imgcollection.map(function(img) { return img.clip(fc); });

  var withCloudiness = combine.map(function(image) {
    var cloud =       ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
    var cloudiness = cloud.reduceRegion({
      reducer: 'mean',
      geometry: combine,
      scale: 30,
    });
    return image.set(cloudiness);
  });

  var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 10));
  print(filteredCollection);


texts:
I have a large set of images that I want to stitch the together and orthorectify. They have been taken from a small aircraft. The camera of the aircraft is tilted, thereby, although gps data is available, the image position and the gps position are different. I have tried some software such as Microsoft ICE and PTgui with no luck.

How could I proceed?


texts:
Does anyone know if there is a plugin to make cartograms for leaflet? I have sets of geotagged Twitter datas and I'd like to distord the map according to the density of tweets of every area of the world.

As cartograms are mostly changing the tiles I'm not sure if Leaflet is the most interesting plugin to do so (as it's mostly done to display datas on a static map layer like an OSM one), but as I begun to use it a lot for other purposes it would be easier to implement in my code.

Another approach would be to find a plugin to transform directly the OSM layer without using Leaflet, but I don't know if that exists too.


texts: How to Load Relational Database Format Data into PostGIS?

texts: I want to load some sample Relational Database Format Data  to POSTGIS. Are there any tools or methods to do this?


texts: Does Calculate Value (Model Only Tool) work correctly in tools run in Batch?

texts:
Am I correct that if my model has a model-only tool (Calculate Value), and then I run the model as a batched process, it will not run the Calculate Value expression?

I wrote a model to upload, define, project (and clip) XY data from excel files. I had thought to use iterate tables, but it wasn't working. (I assume due to having to choose which worksheet within xlxs). So I thought to do a batch, but this seems to leave that annoying $ that arc appends to table names (and is an invalid character in most tools...)

My guess is that the Calculate value that removes my $ is not being run in the batch. Ideas? Confirmations?


texts: ogr2ogr - gpx to PostGIS has wrong timezone adjustment


texts: Trying to export vector shapefile as eps file from ArcMap to CorelDraw?

texts:
I am trying to export a complex shapefile from ArcMap 10.4 into Coreldraw X6. I have already exported other shapefiles by exporting the map as an eps. file and then importing it into CorelDraw which works fine. This one layer won't work though and I get the error 'An error occured while processing .eps. The file cannot be opened.' It is a shapefile with ~12000 small polygons so I have tried removing the tiniest ones and now have ~7000, but this still won't work.

Has anyone experienced this or does anyone have any advice?

I also tried saving it and importing it as an ai. file but that won't work either.


texts:
I'm having an issue with the Layer list popup menu in my WAB apps being slightly transparent, which can make the layer names harder to read.  I would like to set the background opacity to 100% so that the map cannot be seen behind it.
I'm creating the apps w/ WAB Dev Edition and using the foldable theme.  I can't seem to find the setting for the menu opacity anywhere - just transparency settings for the individual layers.
Where can I find and change this setting?  None of the other widgets I'm using seem to have this issue.


texts: Need help for getting a high resolution satellite image for digitising features in ArcGIS?

texts: I want to digitise some places of my city in ArcGIS. I found ArcGIS 'World Imagery' layer which I can use it as a basemap for digitisation. But the problem is it works only while the internet is connected. Is there anyway to save the cache of the layer at a particular zoom level so that I can work on it offline even after I restart the computer. If no, please tell me any other way- may be with google or bing imagery, I need to digitise important features of the city. So i require resolution at house level.