texts: Area-weighted average raster values within each SpatialPolygonsDataFrame polygon (R)

texts:
I'm working in R. I wish to find the average value of all rasters cntrystack overlaid by each polygon within a SpatialPolygonsDataFrame cntryshp. The raster object is actually a RasterBrick. So I want to:

extract all raster cell values, by individual polygons within cntryshp, for EACH layer in the RasterBrick,
average values ACROSS raster cells associated w/ each polygon in cntryshp, for each layer in the RasterBrick, but with the average weighted by the proportion of each raster that falls within the relevant polygon,
bind those data to the cntryshp polygon IDs that identify each polygon, to make a data.frame with 1 row per polygon, and one column per layer of the RasterBrick.

A picture of the polygon over the top layer of my raster brick:

At the moment, I'm successfully doing this, except without the weighting. So my code goes as thus:

myfunc <- function(x) { x.mean <- colMeans(x, na.rm = FALSE, dims = 1) }

r.vals <- extract(cntrystack, cntryshp)

r.mean <- lapply(r.vals, FUN=myfunc)

mean.matrix <- do.call(rbind, r.mean)

countrydata <- data.frame(cntryshp@data$ID,mean.matrix)

However, I wish instead to perform a weighted average, in that function, the weight being the proportion of the raster that falls within each polygon. So to be absolutely clear, consider the bottom right polygon (call it poly1) in my posted graphic above. It falls over a number of rasters: [2,1], [3,1], [4,1], [4,2], [5,1],[5,2]. The raster that holds the largest chunk of poly1 is [4,1]. Then [3,1], which holds slightly less of poly1. [2,1] and [5,2] both hold quite small chunks of poly1. [4,2] and [5,1] both do hold a TINY bit of poly1, but the areas are close to negligible. So I want to average the values of each of these 6 rasters, but weighting them according to proportion under poly1, such that the value of [4,1] is weighted most heavily, and the values of [4,2] and [5,1] are weighted very, very little.

One way to do this would be to use instead the weighted mean function:

myfunc <- function(x) { colWeightedMeans(x, w=wghts, na.rm=FALSE, dims = 1)  }

But that requires extracting the proportion/area of each raster w/in each polygon ("arealist"), I think ideally as a list, so that arealist[[i]] gives the vector of raster areas for rasters that fall within polygon i, ordered in the same way that the rasters values are ordered in r.vals, after extracting. I'm struggling with this (perhaps a raster layer needs to be transformed into a polygon, and intersections taken?), and worried that it might not even be the best approach.

An addendum: I figured out that a built in R functionality does this task perfectly. I had simply overlooked it. So now my code is as follows (actually shorter than the original code):

polymeans <- extract(cntrystack, cntryshp_rdf, weights=TRUE, fun=mean)
countrydata <- data.frame(cntryshp@data$ID,polymeans)


texts: Add GEOJSON layer in tilemill EACCES symlink Error

texts:
I'm trying to add a GEOJSON layer to Tilemill [ Version 0.10.2 (v0.10.1-293-g697c86c) ], installed on ubuntu 14.04.  I've confirmed that the ID used is unique and attempted to load the file into my project via 'Add Layer'.

I've tried to select the file on the local file system, and after I select and click on Save, I get the error, EACCES symlink [filepath].

This is intended to replace an existing .shp labels file.  In Tilemill the .shp file is successfully assigned and displayed, but when I put the newly created .json file into the layers directory where the .shp file is and give the .json file the same permissions I still get the EACCES symlink error.

I've even attempted to serve the file via http on localhost.  I get see that the file is retrieved with a 200, but the same error is seen.

Is this a bug, or am I missing something?

BTW - the geojson I'm trying to load is the output of Dymo, so I'm creating the format myself.  I've confirmed that the target json file can be loaded by python's json library, so it seems there is not a problem with the format.


texts: Is it possible to view an array as label in QGIS?

texts:
I tried to make a GeoJSON file where some properties have arrays:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              4654842.272603876,
              5135720.4450117
            ],
            ...
            [
              4653796.093490953,
              5124881.326846274
            ],
            [
              4653793.42338882,
              5124889.5345974155
            ]
          ]
        ]
      },
      "properties": {
        "id": 72757,
        "parent_id": 69911,
        "common_id": "61",
        "common_name": "81-61",
        "has_children": false,
        "shape_type_id": 4,
        "shape_type_name": "Precinct",
        "value": [
          "54.46",
          "29.34",
          "9.15",
          "4.93",
          "0.47"
        ],
        "color": null,
        "title_location": "Precinct: 81-61",
        "indicator_name_abbrv": [
          "Margvelashvili (Georgian Dream)",
          "Bakradze (UNM)",
          "Burjanadze (United Geo)",
          "Targamadze (CDM)",
          "Liluashvili (Init Grp)"
        ]
      }
    }
  ]
}

However, when I try to use them as labels they don't show up :

So, is it possible to view arrays as labels in QGIS ?


texts: Cannot install any plugins: please see error

texts:
Using QGIS version 3.4.4.
No issues with version 2.18.28

When attempting to install any well credited plugin, I'm presented with the following error.

`Couldn't load plugin 'nominatim' due to an error when calling its classFactory() method

TypeError: 'QVariant' object is not subscriptable
Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python\qgis\utils.py", line 335, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/Users/jamesbrett/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\nominatim\__init__.py", line 27, in classFactory
    return nominatim(iface)
  File "C:/Users/jamesbrett/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\nominatim\nominatim.py", line 54, in __init__
    self.myLocale = locale[0:2]
TypeError: 'QVariant' object is not subscriptable

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

Python Path:
C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python
C:/Users/jamesbrett/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/jamesbrett/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins
C:\Program Files\QGIS 3.4\bin\python37.zip
C:\PROGRA~1\QGIS3~1.4\apps\Python37\DLLs
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib
C:\Program Files\QGIS 3.4\bin
C:\PROGRA~1\QGIS3~1.4\apps\Python37
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\win32
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\Pythonwin
C:/Users/jamesbrett/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/jamesbrett/AppData/Roaming/QGIS/QGIS3/profiles/default`

When closing the window and return to the plugin installer dialogue, the plugin is shown to be installed but the checkbox is unchecked. When checking it, the error returns.

Here is a question with the same issue but I do not know how to resolve the issue.


texts:
I have about 250 rasters (aerial imagery) to merge together in QGIS 3.2. Approx 4.5 gigs in size.

What's the easiest way to merge them? I have tried to "merge" tool a lot but getting various errors such as:

ERROR 1: D:/Desktop/1/1233.tif, band 3: An error occurred while writing a dirty block

ERROR 1: D:/Desktop/1/1233.tif, band 1: An error occurred while writing a dirty block

ERROR 1: JPEGPreEncode:Strip/tile too large for JPEG

ERROR 1: WriteEncodedTile/Strip() failed.

The rasters are in a .tif format.


texts: Why are labels overlapping in ArcGIS Desktop?

texts:
I'm not new to GIS, but just learning ArcGIS 10.

I have been led to believe that the standard labeling engine doesn't allow any overlap: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Positioning_labels_with_the_Standard_Label_Engine/00s800000020000000/

however, I am seeing significant amounts of overlap between labels in just one layer. This is the image I get:

This is within a Data Frame with a fixed scale (1:2500) and no reference scale . I wondered whether that had anything to do with it, but playing around with either of these parameters doesn't seem to have any effect. Is there anything I'm missing in the label manager options - I have checked and double checked that the "Place Overlapping Labels" check box is off. When I increase the buffer size, I get fewer labels, but some of them still overlap.


texts: Standardizing number of jobs and sales per state for choropleth map?

texts:
I want to show the change in jobs and sales per state from 2006 to 2015.

Since I will be making a choropleth map, does this mean that I need to normalize the job and sales numbers before I make the map?

If so, what should I divide them by?

Population?

State area?


texts: I am using the osm2pqsql utility to convert my osm file to PostgreSQL but it always stops showing the following error. How should I solve this problem?


texts: Change "more Info" in ArcGIS Online popup

texts:
I would like to change the text for the "more Info" popup which appears in my ArcGIS Online map next to my URL field. I would like it to say "Show image" or something more relevant than "more info".

Is there any way to change it?


texts: Where can I find base data shapefiles for whole US?

texts:
I am looking for the base layers for whole U.S. - roads, hydrological lines etc with flexibilities of changing symbology and labels.

TIGER is a good source but it is at county level. I am running on a short deadline right now, so want a quick fix. Can anyone suggest good source where I can find one layer for whole country? Or merged TIGER Census shapefiles for the whole US?


texts: Accessing ArcGIS layers from separate thread using ArcGIS Engine?

texts:
I've got an ArcEngine application, referencing a MXD for it's data (and map display).  I'm running a query against all layers, so can take quite some time.

I would normally run the query in a separate thread, which would allow me to return the results as they come in for each layer, rather than wait till all are layers are queried, then return a giant result set.

The limitation however is that I can't access any of the existing ArcObjects in my second thread, due to the STA threading model of ArcObjects.

The only way around it at this stage, appears to be to re-create my data environment from scratch in the second thread, simply for querying.  The problem with this approach is that the layers have joins and attribute queries against them.  For a number of reasons, it's just not a simply exercise to go this route.

I've looked at the IXMLSerialize interface, hoping to find some object, which would allow me to bridge between threads, which I could deserialize/serialize, but nothing sprang out to me.

Any other ideas?

UPDATE:
Because it was so important that I query against the layers, not the datasource, at this point, I save the search layers to .lyr files, then in my second thread, load them using the ILayerFile interface.  This then gives me access to the layers with joins/relates/definition queries intact. The problem at this stage however, is that the ILayerFile.Layer call is very slow.  So am looking for a way to speed that up.  The threaded search does however perform wonderfully, and has improved the UI greatly!


texts: Merging selected features in QGIS: some "borders" remain

texts:
I'm trying to understand (and solve) the following issue in QGIS (2.16.3): I have a shapefile with hexagons which I want to merge. Doing so works, but certain "borders" (not sure about the correct term, sorry) inside the polygon remain, and I cannot figure out how to remove them (see below).

input:

output:


texts: How do I split polygons into line segments?

texts: How can I split my polygon file into "segements"? The lines should be broken up like this picture shows:


texts: GDAL Translate with WMS Source tears raster

texts:
I am using gdal_translate.exe to get a PNG file from a WMS server. Depending on the image size the resulting PNG files shows some sort of fragmentation. The type of fragmentation varies depending on the size and aspect. It seems to happen more with flatter images than square ones but I can't get a consistent handle on the issue.

This happens with the GDAL precompiled binaries from 1.8.0 and also my own local build at 2.1.2 and also the OSGeo4W binaries. There are no problems getting the correct PNG file using a browser, CURL.EXE or QGIS.

So I suspect that something is wrong in GDAL.

Note that -outsize must be set to match the WIDTH and HEIGHT parameters.

THIS IS THE TORN IMAGE (attached)

gdal_translate -outsize 1025 604 -of png "WMS:http://maps.six.nsw.gov.au/arcgis/services/public/NSW_Imagery/MapServer/WMSServer?SERVICE=WMS&REQUEST=GetMap&VERSION=1.3.0&CRS=EPSG:4326&FORMAT=image/png&STYLES=default&LAYERS=BestImageryDates&WIDTH=1025&HEIGHT=604&BBOX=-38.35116600,139.23662900,-26.60347900,159.17524268" "c:\temp\gda_212A.png"

THIS GIVES THE CORRECT IMAGE BUT SMALLER

Hmm - as I don't have enough credits to add more than 2 links I cannot add both the sample image and the second query that returns the good smaller image. So instead I added the image plus instructions on how to get the smaller image.

Replace the following parts of the query as follows...

-outsize 1025 604      >EDIT> -outsize 512 302

&WIDTH=1025&HEIGHT=604 >EDIT> &WIDTH=512&HEIGHT=302

"c:\temp\gda_212A.png" >EDIT> "c:\temp\gda_212B.png"

Does anyone know what the cause is?

Further testing results just now... I changed the output format and request to BMP and TIFF and the result is the same, the image is torn.

PART SOLUTION

I have found what appears to be problem in the GDAL code, but am not sure if it really is. For my two cents worth whilst debugging I found the following problems that seem to be 1.3.0 related:

In function GDALWMSRasterBand::ReadBlocks()

The URL for the block comes from AskMiniDriverForBlock(&url, ix, iy);

The returned block data is based on swapped values for

m_data_window.m_x0

m_data_window.m_y0

m_data_window.m_x1

m_data_window.m_y1

I found that the WMS fetches data in blocks m_default_block_size_x and m_default_block_size_y which are 1024. This triggers the block read and messed up blocks.

I think this could be fixed in the code but can be worked around by using a service description to get the data instead.

Here is something that works now...

gdal_translate -outsize 1161 689 -of png "<GDAL_WMS><Service name=""WMS""><Version>1.3.0</Version><ServerUrl>http://maps.six.nsw.gov.au/arcgis/services/public/NSW_Imagery/MapServer/WMSServer?</ServerUrl><CRS>EPSG:4326</CRS><ImageFormat>image/png</ImageFormat><Transparent>FALSE</Transparent><Layers>BestImageryDates</Layers><Styles>default</Styles><BBoxOrder>yxYX</BBoxOrder></Service><DataWindow><UpperLeftX>139.23662900</UpperLeftX><UpperLeftY>-38.35116600</UpperLeftY><LowerRightX>159.04270529</LowerRightX><LowerRightY>-26.60347900</LowerRightY><SizeX>1161</SizeX><SizeY>689</SizeY></DataWindow><Projection>EPSG:4326</Projection><BandsCount>3</BandsCount><BlockSizeX>1161</BlockSizeX><BlockSizeY>689</BlockSizeY></GDAL_WMS>" "c:\temp\gdal_212_SD1.png"

Justin


texts: Selecting Point Features by Distance and Name

texts:
I have a multi point Townland file. I want to be able to remove duplicates within a radius of 3000 meters from each individual point??

I cannot simply remove duplicates as they are from various locations.

E.g. There may be five points called Castletown within a 4000m radius and I only want one.

However there may be another five points called Castletown 40,000m away (As more than one area is called Castletown)  and I also only want to keep one of these within a 4000m radius of itself.

Therefore I would end up with two points called Castletown from two separate locations.

Any ideas whats the best way to do this?


texts: Can I export QGIS tool results to anything other than a shapefile or DBF?

texts: I am trying to use the QGIS 'Split Lines at Points' tool, since my organization does not have the package to unlock the ESRI version of the tool. I am able to access my data in an enterprise ESRI Geodatabase just fine, however the only outputs listed for that tool are shapefile and dbf. We use feature classes in our database organization, and my supervisor is wary of converting the data from a shapefile back to a feature class. Is there any plugin I can add to QGIS to allow the tools to output the data to feature classes?


texts:
I am trying to publish a draft service definition using arcpy.mapping.CreateMapSDDraft and am getting the referenced Assertion Error. I've tried the responses to users with similar questions with no luck.
Here is the code I am working with:
# Import system modules
import arcpy, os, sys
from arcpy import env

# Set local variables
wrkspc = r"\\intruder11\grc\sw\core\dev\gp\hosted_services_sp"
mapDoc = arcpy.mapping.MapDocument(wrkspc + 'NoMap.mxd')
outputFolder = r"\\intruder11\grc\sw\core\dev\gp\hosted_services_sp"
sddraft = mapDoc + service + '.sddraft'
summary = 'Test of publishing automation'
tags = 'CNR Southwest, dev, test, acc'

# Sign-in to Portal for ArcGIS
logger.info('Signing into built-in ACC account')
arcpy.SignInToPortal_server("--","--", "https://--")

### create service definition draft
analysis = arcpy.mapping.CreateMapSDDraft(mapDoc, sddraft, service, 'MY_HOSTED_SERVICES',con, True, None, summary, tags)

### Sign out from ArcGIS Online
arcpy.SignOutFromPortal_server()


texts: Export multiple bands from netCDF into GeoTIFFs

texts:
I am working with a NetCDF file that has multiple subsets and each subset has multiple bands. I would like to take one subset and export each of the 28 bands into its own GeoTIFF. I thought that gdal_translate could do this via the -sds flag, which is described in GDAL documentation as "Copy all subdatasets of this file to individual output files. Use with formats like HDF that have subdatasets", however I guess this doesn't carry over to bands as well. This command:

gdal_translate -of GTiff -s_srs epsg:4326 -a_nodata = 0.0 -sds NETCDF:"oscar_vel2020_180.nc”:u temp.tif

generates the terminal prompt >.

How might I export each of the bands from the subset of the NetCDF file into its own GeoTIFF?


texts:
Already read Merge Multi-Band Geotiff with different datatypes (with gdal) but can't apply that solution so please do not flag as duplicate.

I have near 3.000 geotiff files that I need to merge into just one.

this is what I done so far:

gdalbuildvrt -input_file_list list_of_files.txt -overwrite -addalpha images.vrt

and then

gdal_translate -tr 0.000170 0.000170 -r cubicspline -of GTiff images.vrt test.tif

I've got this warning:

Warning 1: Unable to export GeoTIFF file with different datatypes per
different bands. All bands should have the same types in TIFF.

I can get rid of the warning when remove -addalpha from gdalbuildvrt.

I'm a complete newbie on this so I ask: What is the implication of the -addalpha parameter and that warning if I choose to keep it?


texts:
How can I highlight nearest(certain distance) unique duplicates? Actually I tried a process in ArcGIS but only duplicates are pulled out from the entire file without considering distance. How can I do the process in QGIS.

.


texts:
I have a local system coordinates (Pi) of a satellite and I want to calculate observer latitude and longitude from the local coordinates of the satellite. I did this operation manually, I didn't use any software.
The following image shows the local system and the global system.

n(unit normal vector) is defined:

For this purpose I did these steps:

I found the zenith angle (alpha) of the satellite in the local system.

I applied the transformation matrix to Pi using the zenith angle. Transformation is counter-clockwise direction. I used CCW direction as the positive angle. The result vector Pi_new.

I thought the result vector Pi_new (finding in 2. entry) orientation is the same as normal vector (n) orientation. So to find the unit normal vector, I divided the
Pi-new / |lenght_Pi_new| (normal unit vector is Pi_normal)

4.To find the latitude and longitude of the observer I use normal vector (n) and Pi_normal.I thought they should be equal.

longitude = arctan(Pi_normal_y,Pi_normal_x)
latitude = arcsin(Pi_normal_z)

I did all these steps but, when I tested I did not get the true latitude and longitude values.
Where am I doing wrong?
*I get the local system as the left-hand system.
Edit:
The output coordinates should be
lat: 40.0023 and lon: 32.9024

but my calculation result
lat: 27.94597928  and lon :1.31664029

There is huge differences.
Also my input (satellite in local system) is
[[ 24653272.96068735]
 [-13050866.8145796 ]
 [ -1065587.51414906]]


texts: How to get elem_info array from WKT or JTS Geometry in Java?

texts:
I have Well Known Text (WKT) string with geometry.
Now I want to create Oracle Mapviewer SdoGeometry, which needs sdo_elem_info array as one of arguments which is an array of doubles.

How can I extract elem_info array from WKT String?
I've read that geometry into Geometry object from Vivid Solutions JTS library, but I can't find a way to get elem_info from that either.

Anyone can help?


texts:
I am involved on a project which is about Generating 2 maps via running macros in Arcmap 9.3.

It works fine when user push the "run macro" button but now I want to call arcmap from a web page.

It opens arcmap file successfully but I was wondering is there a way to run macro without clicking? This project is going to be part of web application and must work automatically.

Thank you


texts: Use GeoServer as source of vector layer in OpenLayers 3

texts:
I am using PostGIS and GeoServer and OpenLayers 3 (OL3). I connected PostGIS table to GeoServer and published my layer.  I am using following code to get the layer as vector layer:

var vector = new ol.layer.Vector({
    source: new ol.source.GeoJSON({
        projection: 'EPSG:4326',
        url: 'http://localhost:8000/geoserver/cite/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=cite:myfeaturesTbl&maxFeatures=50&outputFormat=application/json&'
    })
});

I need to read my feature as vector layer because I want to select them and then get information of each feature or move them on the map. When I use this code it is not show anything.

Is there any problem in my code?

Can you help me get a vector layer in OL3 with a GeoServer WFS?


texts: Reclassifying features before exporting data driven page

texts:
I need to find a way to reclassify/refresh graduated color symbology in between exporting data driven pages.

Background: I have spatially joined Block Groups (BGs) to 1000 facilities (and facility ID) within 10 miles. There are multiple instances of the same BGs because they are located in close proximity to multiple facilities.

I prefer to use the standard deviation classification method. Currently, when Data Driven Pages are enabled and I classify the data - the value ranges are applied are relevant to what is displayed (this is good). The problem arises when I toggle to the next page and this is not re-calculated, but rather continues to use the value ranges for the previously displayed BGs.I cannot display all at once, because the multiple occurrences overlap. Also, for one particular facility a BG might represent a lower value while for another the same BG might represent a larger value.

I export data driven pages in very basic python script. One answer would be to insert code that applies symbology from a different layer before exporting the page. This process works manually in Arcmap: Apply Symbology > refresh > export page...etc. I have limited experience with python, and the 2nd to last line I have added below isn't making it happen (but works in a separate script).

This is the python I am currently working with:

import arcpy
mxd = arcpy.mapping.MapDocument(r"E:\Templates_MXD\BU_5Mile.mxd")
export_path = r"C:\Users\sder\Desktop\Test\\"
field_name = "ID"
df = arcpy.mapping.ListDataFrames(mxd, "Business Analyst Map")[0]
updateLayer = arcpy.mapping.ListLayers(mxd, "BKtoBG", df)[0]
sourceLayer = arcpy.mapping.Layer(r"E:\Templates_Layer\BKtoBG.lyr")

for i in range(1, mxd.dataDrivenPages.pageCount + 1):
    mxd.dataDrivenPages.currentPageID = i
    row = mxd.dataDrivenPages.pageRow
    print row.getValue(field_name)
    arcpy.mapping.ExportToPDF(mxd, export_path + row.getValue(field_name) + ".pdf", resolution=200, image_quality="BETTER")
    arcpy.mapping.UpdateLayer(df, updateLayer, sourceLayer,True)
del mxd


texts: Is it possible to use a raster::click event within an interactive leaflet map?

texts:
Data Download: https://www.dropbox.com/s/f849ra8l6umo9d8/crops.tif?dl=0

I have various raster images of the world that I would ultimately like to plot within an interactive leaflet map; an example is linked above. With basic R plotting, I am able to plot the .tif image and use raster::click to be able to click on any location on the map and have the associated legend value printed in the R console. The code is very straightforward:

library(raster)
library(leaflet)

setwd("C:/Users/administrator.me/Desktop")

crop <- raster("C:/Users/administrator.me/Desktop/crops/crops.tif")

plot(crop)
click(crop)

Again, the result is a raster image of the world that can be clicked and spits out the associated value (like using the identify tool in QGIS):

This is pretty simple, but I would ultimately like to have this raster image as a background layer in an interactive Leaflet map for R. I've plotted it with the additional code:

library(leaflet)

pal <- colorNumeric("Spectral", values(crop),
                    na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(crop, colors = pal, opacity = 1) %>%
  addLegend(pal = pal, values = values(crop),
            opacity = 1,
            title = "Crop Calendar")

This essentially plots the same image as above, just over a basemap and with panning/zooming options. I would like to be able to use the same raster::click event within the interactive leaflet map. I've read a little bit about click events in Leaflet (example here], but it seems like these are mostly for predefined markers or polygons. Is it possible to do this with a gridded raster image? Can I create a click event that gives a pop-up window of the value associated with the grid that is clicked on?


texts: Extracting feature attributes/geometry using PyQGIS

texts:
I am trying to extract co-ordinate and height data from a contour layer imported from a .shp file. I have loaded the file into QGIS and am attempting to access with Python the height values of each feature - Closest vertex Z or Interpolated Z in the picture below.

I am eventually aiming for something similar to what is achieved in this answer https://gis.stackexchange.com/a/41922/172876, but I can't find a way to adapt their solution to work for me. I have looked through several forum posts, but many rely on using rasters to extract z values, or are using points not contours.
I have been able to select features in the contour layer, and access their attributes using:
layer = iface.activeLayer()
for feature in layer.selectedFeatures():
    attrs = feature.attributes()
    print(attrs)

But this only accesses the attributes pictured below, and none of the attributes under the (Derived) dropdown in the picture above, which is where the Z values are.

How do I access and extract these using PyQGIS?


texts: How to clean such a bad shapefile for routing?

texts:
I am trying to clean a shapefile to make it usable for routing with pgrouting, postgis, postgresql.

Here an intersection, it seems to be connected but when we zoom …. It is not:

Same here :

Here (it is zoomed at max) there is 2 ways when it seems to be one, and only one is connected :

I am using assign_vertex_id with different tolerances, and QGIS and PostGIS, but I don't know what to use in it.


texts: I want to edit corresponding postgresql data of a vector from openlayers- how to do it?

texts:
I want to edit the corresponding text/numeric data type columns of a vector through openlayres, when a user click/mouse over a vector. like,

http://dev4.mapgears.com/bdga/bdgaWFS-T.html#

any tutorial page for how to do it? I gone through it , it was difficult for me to understand from its js codes.

-- asked in stackoverflow, but no answers. based on a comment, i posted the same question here.


texts: Diagnosis page Error:Trying to visualize a point data using createviz in cartodb.js

texts:
I had setup the cartodb locally in my machine following this tutorial https://github.com/danseely/cartodb-install/blob/master/DEV-INSTALLATION.md

Everything is installed perfectly: I can visualize a data in my local account, but when i try to import a dataset it's showing an error.

So i created dataset with point data. When i try to visulaize it by using the createviz function in cartodb.js i can only see the basemap but my point data is not visible.

When i go to the diagnosis page (development.localhost.lan:3000/diagnosis)

I'm getting this error

CartoDB extension: error

    Installed version: 0.11.1
    Current supported version: 0.9.Latest version: 0.9.3

Windshaft: error

    Service url: http://localhost.lan:8181
    Full config: {"filter"=>"mapnik", "internal"=>{"protocol"=>"http", "domain"=>"localhost.lan", "port"=>"8181", "host"=>"127.0.0.1", "verifycert"=>false}, "private"=>{"protocol"=>"http", "domain"=>"localhost.lan", "port"=>"8181", "verifycert"=>false}, "public"=>{"protocol"=>"http", "domain"=>"localhost.lan", "port"=>"8181", "verifycert"=>false}}
    error fetching info: A JSON text must at least contain two octets!

SQL API: error

    Service url: http://localhost.lan:8080
    Full config: {"private"=>{"protocol"=>"http", "domain"=>"localhost.lan", "endpoint"=>"/api/v1/sql", "port"=>8080}, "public"=>{"protocol"=>"http", "domain"=>"localhost.lan", "endpoint"=>"/api/v2/sql", "port"=>8080}}
    error fetching info: A JSON text must at least contain two octet


texts: File locked after removing the layer , in Python code for a plugin for QGIS

texts:
I try to remove the layer, to re-generate it with the same name, but when the plug is started the layer is loaded, removeMapLayer does not unlock the file, and it is locked.

I have seen several messages from people who have had the same problem but no solution is currently serving me.

QgsProject.instance().removeMapLayer(layer_.id())
os.remove(url_layer)


texts: repeat skobbler maps horizontally with OpenLayers?

texts:
I'm trying to get skobbler maps to repeat horizontally
if you load this page
http://developer.skobbler.com/demo/skobbler-tiles-openlayers.html
and zoom out, the map will not repeat horizontally

However if you use this link and zoom out
http://developer.skobbler.com/web-tiles/embedded
it works as I want (i.e. it repeats across the screen)

However, I'm using OpenLayers and not an Iframe

Is what I want to do possible with OpenLayers?

Thanks
Mark


texts: How do I extract points from a LineString or Polygon data frame and make a dictionary of point data frames?

texts: I am building some pre-processing routines for 3D visualization of GIS data and I have several line and polygon shapefiles with multiple features. I want to extract the points for each of the features along with the feature attributes and create a dictionary of point geo-dataframes. This I suspect is easy to do, but I am not quite sure how to index into the linestring within a geodataframe.


texts: QGIS Save as Spatialite layer from FileGDB does not add OGC_FID

texts: I have previously (yesterday) been able to save a Spatialite layer from a FileGDB feature class using the Right click/Save as... method in QGIS. Today, however, with no changes to any software (tried on 2 machines), the creation of the Spatialite layer today does not appear to include the OGC_FID field.


texts: create hub distance - restricted to state boundaries

texts: I am trying to calculate hub distance for enumeration areas within a shapefile of a country using QGIS 2.18. The create hub distance in QGIS calculates this based on the distance to the nearest hub. As a result some hub lines cross over into another state from their original state. Is there any way I can do this without having having to create a shapefile for each of the states as they are more than 40 states in my shapefile.


texts: Using ArcGIS field calulator with VB Script

texts:
My code looks like the following:

if [LWS50yrp] <136.5 then
     [DMG_50y] = 0.1
elseif [LWS50yrp] >=136.5<154 then
     [DMG_50y] = 0.15
elseif [LWS50yrp] >=154<166 then

...

elseif [LWS50yrp] >=400<500 then
     [DMG_50y] = 1.8
elseif [LWS50yrp] >500 then
     [DMG_50y] = 2
end if

Calculation runs all fine. However all values from [LWS50yrp] exceeding 136.5 are assigned to a value of 0.15!

I can't find the mistake, must be a small twist, but I'm only a beginner and have no clue.


texts: Adding attributes to a feature using ChangeAttributeValues in Pyqgis

texts:
I have a workareas polygon layer with X number of features and a buffer CSV layer with Y number of entries which looks like the table below (one id field and a couple other fields).

id1  | A1 | B1 | C1 | D1 | E1
id2  | A2 | B2 | C2 | D2 | E2
id3  | A3 | B3 | C3 | D3 | E3
id4  | A4 | B4 | C4 | D4 | E4
id5  | A5 | B5 | C5 | D5 | E5

What I want to achieve is to duplicate every workarea feature Y times and copy the CSV table attributes. I've pasted my code below, but the changeAttributeValues does not work for me.

workareas_mem = QgsVectorLayer("Polygon?crs="+workarea_lyr.sourceCrs().authid(), 'polygon', 'memory')
prov = workareas_mem.dataProvider()
prov.addAttributes([field for field in workarea_lyr.fields()])
prov.addAttributes([field for field in buffer_tbl.fields()])
fld_names = [field.name() for field in buffer_tbl.fields()]
workareas_mem.updateFields()

if True:
    for i, workarea in enumerate(workarea_lyr.getFeatures()):
        for buffer in buffer_tbl.getFeatures():
            new_feat = workarea
            prov.addFeatures([new_feat])
            for fld in fld_names:
                prov.changeAttributeValues({new_feat.id():{prov.fieldNameMap()[fld]:buffer[fld]}})
        feedback.setProgress(int(i* total))

The code successfully creates the duplicate workarea features and I end up with X*Y features in total. But the attributes from the CSV table does not write properly. Attached below is a sample output from this code.


texts: Splitting shapefile into separate shapefiles using ArcGIS Desktop

texts:
I have a multipart shapefile, which I would like to split into many individual shapefiles (by FID for example).

Is there some way to do this in ArcGIS 10.x without using python?


texts: QGIS - how to move layer files while project is open

texts: I want to move layer files (within browser directory structure) while my QGIS project is still open - therefore allowing the project to update its relative link to that file. ARCGIS does this, but not QGIS? Is there a workaround?


texts: Assigned projection and an unprojected WGS84 geographic coordinate system

texts:
I am to be submitting a shapefile for review but the first two requirements have me confused. Are these requirements conflicting or am i just not understanding the question?

Requirements:

The shapefile must have an assigned projection w/accompanying .prj file.
The shapefile must use an unprojected (geographic) WGS84 geographic coordinate system.

I will be converting the original data from a .kml file which is automatically set to the WGS 1984 coordinate system when exported to the shapefile. Is there a way to assign a projection?


texts:
I am making a plugin. I have two QML files for the same layer and column. I want the plugin to combine the two existing files into a new QML file and then I will select which one symbol for each value I want.

The problem is, I don't know how to combine two QML files into new one. Any help?


texts: Gdal: how to get the max and min altitudes of my topographic raster?

texts:
Given a .tiff topographic raster image.

How to get the altitudes of the highest and lowest point/pixel ?


texts:
Inserting attribute table into QGIS Composer is limited to 99 rows?
Or am doing something wrong?

Stems from
How to import or incorporate formated Excel tab into QGIS print composer


texts: Find the vertices of an ROI in XROI using IDL

texts:
I need to find the vertices of an ROI shape created in IDL's XROI, so that I can redraw the ROI on multiple images. I know you can find the centroid, area and perimeter with the GEOMETRY command, but I can't find anything about vertices.

Is there a way of finding the coordinates of an ROI's vertices in IDL?


texts: How to simplify a polygon down to X amount of points?

texts: If I have a polygon that has 100000 points, how can I calculate what the tolerance should be to simplify it down to 1000 points using ST_SimplifyPreserveTopology?


texts: Adding scale bar to raster overlayed on ggmap Google Maps object

texts:
I am trying to plot a predicted surface (raster layer) over a ggmap object that contains base imagery from Google Earth. I have produced the map that I want however, I am struggling to add the legend scalebar.

To add the raster over top the ggmap object, I used the inset_raster function and I am not sure where to map the aesthetics to produce the proper legend based on the color scale I assigned the raster.

Below is an example of my code. Does anyone have advice on how to add the legend based on the inset_raster object?

rsfmap <- ggmap(get_googlemap(center = c(lon = x, lat = y),
                              zoom = 9,
                              maptype = "hybrid"),
                              extent = "device", darken = 0.05, legend='right')

rsfmap+
  geom_point(aes(x=coords.x1, y=coords.x2), size=1, alpha=0.75, data=mydata) +
  inset_raster(as.raster(ggras, col=(topo.colors(500, alpha=0.4))), xmin =
               ggras@extent[1], xmax = ggras@extent[2],ymin = ggras@extent[3],
               ymax =  ggras@extent[4], interpolate = TRUE)


texts: How to round pixel values of a raster in QGIS?

texts: I have a raster with float number values for example: 2.544459104537964 and I want to round these up to integers. I have seen other posts on ArcGIS like this one but when I use Raster Calculator in QGIS, Int(myraster + 0.5) only makes the first value (which is 0) to 0,5 and the last one (which is 100) to 100,5. The other pixels stay float.


texts: Automating Define Projection within Python script tool?

texts:
How would I automate the Define Projection tool so that I can just enter my layer each time?

I have it as a script in toolbox, but I am unsure how to set it so I can just input my layer as GetParameterAsText

Here is the script. I want to replace the "forest.shp" with a layer of my own choosing each time

 # set local variables
 inDataset = "forest.shp"
 coordinateSystem ="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"

arcpy.DefineProjection_management(inDataset, coordinateSystem)


texts: Extracting data from ArcGIS REST endpoint?

texts:
I am trying to export the data in geoJSON or JSON format from this REST service I keep getting error request can someone help me with the input querys

http://maps3.arcgisonline.com/ArcGIS/rest/services/A-16/FEMA_100-Year_Flood_Zones_in_the_USA/MapServer/1

I need the geometry and attributes returned

I am using this page to try and query for all of the data

http://maps3.arcgisonline.com/ArcGIS/rest/services/A-16/FEMA_100-Year_Flood_Zones_in_the_USA/MapServer/1/query

I choose the format as JSON and in the Return Fields (Comma Separated): I enter "*"

What else do I need to enter into the QUERY to return this layer?

I entered objectid is not null into the where clause and it returned the cryptic ESRI rings format

http://maps3.arcgisonline.com/ArcGIS/rest/services/A-16/FEMA_100-Year_Flood_Zones_in_the_USA/MapServer/1/query?text=&geometry=&geometryType=esriGeometryPoint&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&objectIds=&where=objectid+is+not+null&time=&returnIdsOnly=false&returnGeometry=true&maxAllowableOffset=&outSR=&outFields=%22*%22&f=pjson


texts: I need to obtain a grid of bare ground elevation (removing all edifices, trees, roads and man-made elements) for quantitative geomorphological analysis from a high resolution (1 m) DSM obtained by orthoimages (not from LIDAR). Is there any suggestion you can provide me to filter my original grid?


texts: ArcMap - Strange Behavior When Trying To Wire IActiveViewEvents in IDockableWindow

texts:
I can't seem to wire up IActiveViewEvents for an IDockableWindow.  I am using a straight out-of-the-box IDockableWindow which inherits BaseCommand.  The snippet below shows what I have added to the autogenerated code.  Problem is, when I try to run ArcMap, it (the main ArcMap GUI) never loads.  Please see attached screenshot to see what I mean.

 public override void OnCreate(object hook)
    {
        if (hook != null)
            m_application = hook as IApplication;

        if (m_application != null)
        {
            SetupDockableWindow();
            base.m_enabled = m_dockableWindow != null;
        }
        else
            base.m_enabled = false;

        IMxDocument mxDoc = (IMxDocument)m_application.Document;
        IActiveViewEvents_Event activeViewEvents = mxDoc.ActiveView as IActiveViewEvents_Event;
        activeViewEvents.ViewRefreshed += new IActiveViewEvents_ViewRefreshedEventHandler(activeViewEvents_ViewRefreshed);
    }

    void activeViewEvents_ViewRefreshed(IActiveView view, esriViewDrawPhase phase, object Data, ESRI.ArcGIS.Geometry.IEnvelope envelope)
    {
        MessageBox.Show("ViewRefreshed");
    }


texts: Generate a single shapefile from various intersection agreements

texts:
I have five vector files wherein I want to determine the spatial agreements among them using intersect. Here is my illustration. (Actually these are forest shapefiles from five different forest maps.)

 A = vector_file_1.shp
 B = vector_file_2.shp
 C = vector_file_3.shp
 D = vector_file_4.shp
 E = vector_file_5.shp

So that I will have a set of spatial agreement from 2/5, 3/5, 4/5 to 5/5.
For example:

1)   2/5 = A∩B, A∩C, A∩D, A∩E, B∩C, B∩D, B∩E, C∩D, C∩E
2)   3/5 = A∩B∩C, A∩B∩D, A∩B∩E, A∩C∩D, A∩D∩E, A∩C∩E, B∩C∩E, B∩C∩D, B∩D∩E, C∩D∩E
3)   4/5 = A∩B∩C∩D, A∩B∩C∩E, B∩C∩D∩E, A∩C∩D∩E
4)   5/5 = A∩B∩C∩D∩E

My question is what geoprocessing technique will I use to generate a single shapefile from each agreement above (except item # 4), supposing I've done the intersection combinations for agreement items 1-3? Will I use intersection again, union or other technique?

(My professor told me to just add a field and put a number code on that field, for each shapefile, unfortunately I forgot the next process he mentioned.)


texts: Summing time spent on effort from point GPS data tracks in grid cells?

texts:
I have one layer of point data representing GPS tracks collected during vessel surveys (an individual survey is indicated as the time between departing and returning to shore and there can be more than one on a given date i.e. when we returned to shore and then went out again). I also have layers of 1 and 2 kmsq polygon grids. I would like to calculate the time spent 'on survey effort' in each cell to adjust for effort in other analyses. I intend to adjust the weighting of each cell for these other calculations depending on the amount of effort spent in the cell, i.e. adjusted density of whale sightings based on the search effort in that cell.

So far I have used a Spatial Join in ArcGIS to give the effort points that fall into a particular grid cell the unique identifier of the grid cell. I then exported the table to Excel to calculate the time spent on effort in each cell using the unique identifier of that cell. This was calculated by summing the times (on a given date and survey) between the points collected along that track in that cell and then using a Pivot Table to sum the total survey effort in each cell by the cell's ID and the survey number.  Unfortunately, I realized after doing this that the total time on effort is different at the two different spatial scales of 1 and 2 kmsq. I believe this is because the calculation is missing the times the survey track crossed into a new cell i.e. the points fall in different cells so the time between these points  was excluded.

Does anybody have a good idea of how I could calculate the sum of these survey effort times without losing the times during crossover from one cell to another?

I think I need to make the point data into line data and extract the times when the line crosses into a new grid cell but I'm not sure the best way to do this.


texts: Exposing ArcInfo level liscenses to ArcView desktops via ArcGIS Server and custom python script tool

texts:
I have a script that uses an ArcInfo level tool.  The script works great on my machine but uses an ArcInfo level tool which the ArcView users don't have.

Is exposing these ArcInfo level tools possible throught ArcGIS server?  What I am having trouble conceptualizing is how my python script tool will call that gp service.  Please help me think through the liscensing and calling the gp service from a custom python script tool on the toolbox.

I have ArcInfo, I added the tool to arcmap's TOC saved the mxd and opened the tool and ran it once.  I then published the tool throught ArcGIS Server.  I am stumbling on how to proceed.  Thanks.


texts: Where to download 10m (1/3 arc-second) DEMs?

texts:
I'm looking for some of the 10m DEM to fill in the white Quad areas. Do you know where of can I find them?  These 10 DEMs we have are actually coming from the NRCS Geospatial Dataway. I believe they were not available yet. Are there some web sites that allow me to download some of the DEM ?

What I want is to find them in 10m. I have check out the USGS National Viewer Map and they have 1,2,3 arc seconds.

Place = Mount Baldy (aka Baldy Mountain), which is located on Cimarron Range of the Sangre de Cristo mountains of southern U.S.'s New Mexico state.


texts: Does calculating statistics in ERDAS change raster values?

texts:
I've been under the impression that calculating or recalculating statistics in ERDAS doesn't actually modify an image, but I notice when I do this that when I look at the image file in windows explorer, it says the image was modified when I calculated the statistics.  Does calculating statistics actually change the image values in any way?

The files I've been using have been mostly tiffs, but I think it does the same for .imgs too.


texts:
I am using QGIS: Raster/Extraction/Clipper

I have a raster image (Sentinel 2 satellite 3 bands combination, total area: 1,000,000 hectares, file size 740 MB) from which I clip a certain area, which accounts to 1% of the total image.

When I save it, the size of the clipped raster is 700 MB (only 1,000 hectares),

How can I reduce the clipped file size?


texts: How to get pixel size/cell size of raster layer?

texts:
After much searching, I've still not been able to figure out how to perform what seems like a simple function: I'd like to get the pixel size/cell size of a raster layer using PyQgis. Is there any native way to do this through QgsRasterLayer or a helper class? Or will I need to read the file in using GDAL functions, e.g.:

geotransform = dataset.GetGeoTransform()
if not geotransform is None:
    print 'Pixel Size = (',geotransform[1], ',',geotransform[5],')'

Thanks for your help.

Alex


texts: Basic Python Reclass question in Field Calculator

texts:
I'm trying to reclass values in field calculator using Python.

I want any value >= '5' to remain the same (return same value). Anything else to be = '0'.

Yes, DIST_EDGE are numbers.

I'm getting an invalid syntax error in Line 1 (see below)

Here's what I came up with (new script following comments):

Error Message:


texts: Editing snapped polygons in QGIS: is releasing of snapped vertices possible?

texts:
I have made a layer with polygons that are snapped to each other. Now I would like to edit this layer and create some extra plots in between them.

Is it possible to "unsnap" the vertices? I am using QGIS 2.14.11.


texts: How does majority algorithm works in ArcGIS

texts:
I'm trying to understand how does the majority resampling works in ArcGIS, they write that

"The MAJORITY option, which performs a majority algorithm, determines the new value of the cell based on the most popular values within the filter window. It is mainly used with discrete data just as the nearest neighbor method; Majority tends to give a smoother result than Nearest."

I am talking about the Majority Filter geoprocessing tool from Spatial Analyst extension.

What if there are two most popular values?

Is there a way to resample a raster if for example 50% is covered by a value?


texts: Mercator projection: formula for latitude based on pixels?

texts:
I am basing this question on knowledge from Mercator projection: problem with latitude formula

In that question/answer, the formula is given for working out how many pixels above/below the equator a given latitude will be, assuming you know the total size of the map. The formula is:

Pixel Shift = Width / (2 * pi) * ln(tan(pi / 4 + (Latitude / 2) * pi / 180))

For example, on a 5000 pixels wide map, the latitude 38 would be 571.36 pixels above the equator. I have built a script using this formula and it works perfectly.

What I would like is a similar formula, but which provides a latitude based on a given pixel shift. So I could feed it the answer from above (571.36) and it would give me the answer 38.

I have spent some hours now trying to figure out (based on scraps of remembered high school math) how to balance and switch the formula around but it's beyond me.

Can anyone help me with this?

Here is the code I am trying, based on Russell's answer:

function getPixelShift ( pxWidth, map_latitude )
{
    var pi = Math.PI;
    var pixelShift = (pxWidth / ( 2 * pi) ) * Math.log( Math.tan( pi / 4 + (map_latitude/2) * pi / 180 ) );
    return pixelShift.toFixed(2);
}

getPixelShift (5000, 38); // Returns 571.36

function getLatitude ( pxWidth, pixelShift )
{
    var pi = Math.PI;
    var a = 1 / ( pxWidth / (2 * pi ) );
    a = Math.exp( ( pixelShift * 2) * a );
    var lat = (Math.asin( ( a - 1 ) / (a + 1) ) ) * ( pi / 180.0 );
    return lat;
}

getLatitude (5000, 571.36); // Returns 0.01157551908131

I'm probably implementing what you suggested wrong, sorry! I'm needing the result of getLatitude to be 38.

What am I doing wrong?


texts: Need to get distance matrix for large set of points on-the-fly

texts:
I do not know much about GIS and I do not know much of the jargon. Hopefully you can help me a bit here because I am really unaware of the services and software available.

My problem is as follows: I need to get a distance matrix for a large set of points (up to ~1000 points) all within one city. I need the driving distance (the driving time would be useful as well) between all pairs of such points (so can be up to a million pairs.)
Moreover, I need to do this within seconds.

I am not sure if this is actually possible, if it is not then I'll have to reduce the set of points.
I do know what the city is in advance, but I do not know in advance what the set of points are. Once the set of points are input into the system, I want to be able to get the distance matrix within seconds (~1-2 seconds).

Some ideas I had so far:

Use an API. I did find Google's Distance Matrix API and MapQuest directions API, but neither can handle over 25 locations.

Maybe there are APIs which can handle this?

Store the road network ahead of time and calculate the distance matrix locally. This seems like a more feasible solution.

But where do I get that data and what software/libraries can I use to generate the distance matrix from that data?


texts: Single Time Slider for Multiple Data Frames--ArcMap Desktop 10.3

texts:
This question has been asked before, but never really answered: Is it possible to have a single time slider control that animates over multiple data frames? Reason being, I have one main map and 4 inset maps with the same time-aware dataset (see image with main frame activated). I can only export the animation with one frame following time. I would like the end result to be sequential images.  Has ESRI evolved some way to control multiple frames, or can you suggest a method that would get the same result?

I am aware of the questions at these stackexchange links (How to apply ArcGIS time slider on Multiple Data Frames at once? [3-years old and no solution] and (ArcGis using the time slider to animate multiple data frames [3-months old and no answer]). However, it seems there is a demand for an workaround here, so I would like to keep this question on the boards, though it means reposting the question.


texts: Check if Google Earth Engine Python API is initialized

texts:
In Google Earth Engine Python API, you need to start with an initialization :
import ee

ee.Initialize()

Is there a way to check if the ee has already been Initialized ?


texts: Why would labels not display in dynamic arcserver service?

texts: I’ve got a dynamic map service running ArcGIS Server 9.3.1 that will periodically stop displaying labels for some features.  When you access the service either through a web application, or ArcCatalog, the features draw, but the labels for the features do not draw.  And it’s not consistently the same features.  One time it could be the addresses, the next time it could be the streets.  Has anyone seen this behavior before, and if so are there any recommendations on ways to fix this behavior?


texts: Using remote controlled jeep to gather data points?

texts:
Is it possible to use a small remote controlled jeep or truck - 18" long or so - with some kind of GPS device attached to it to gather data for a detailed topo map of a 2 acre parcel that I just bought?

I was thinking of putting the GPS device in the car and driving it all over the property.  There is a lot of mesquite on the property, weeds, etc.


texts:
I am relatively new to ArcGIS and I had a question. I am working with a data set relating to cities and towns.
I am trying to find cities that have the same name as towns. Currently, the syntax I am using is such: "CSD_NAME" LIKE "NAME". But this is not working, what am I doing wrong?


texts: Joining raster layer and vector (points) in ArcGIS 10

texts: I have a raster layer (land cover) and a vector point layer, and I want to calculate the value of each point in terms of the type of land cover so that it is added to the attribute table. I tried to convert the raster to vector so that I could use the intersect tool or perform a spatial join, but the raster layer is too large (over 2MB) so the conversion failed. I also tried splitting the raster so that it was broken into smaller tiles, but the processing time was taking forever. Does anyone know an alternative way to do this?


texts: Display GeoJSON data on map using django-leaflet?

texts:
I want to display GeoJSON data on map in geodjango.

When I use the below code only GeoJSON data like:

{
  "type": "FeatureCollection",
  "features": [{
        "type": "Feature",
        "properties": {
          "model": "pft.pft"
        },
        "id": 1,
        "geometry": {
          "type": "Point",
          "coordinates": [76.996944, 31.58]
        }
      }, {
        "type": "Feature",
        "properties": {
          "model": "pft.pft"
        },
        "id": 2,
        "geometry": {
          "type": "Point",
          "coordinates": [76.670778, 31.487083]
        }
      }, {
        "type": "Feature",
        "properties": {
          "model": "pft.pft"
        },
        "id": 3,
        "geometry": {
          "type": "Point",
          "coordinates": [77.253, 30.651722]
        }
      },

is getting displayed on my webpage, but I want to display it on map.

template/map.html

{% load leaflet_tags %}

    {% leaflet_js %}
    {% leaflet_css %}

    {% leaflet_map "main" callback="main_map_init" %}

Plant Functional Types

<button type="button" onclick=>Actual plant functional type</button>
<button type="button" >Potential plant functional type</button>
<button type="button" > Plant functional types likely to change</button>
<button type="button" > Upload species data</button>

views.py

def geojson_provinces(request):
    conn = psycopg2.connect(dbname="geodjango",host='localhost',user='postgres', password='postgres', port=5433)
    dict_cur = conn.cursor(cursor_factory=psycopg2.extras.DictCursor)
    cur=conn.cursor()
    res = dict_cur.execute("SELECT ST_AsGeoJson(geom) AS json FROM pft_pft LIMIT 10;")

    points=dict_cur.fetchall()
    print (points)
    json_res  = []
    for row in points:
        json_obj  = row['json']
        print (json_obj)
        json_res.append(json_obj)
        print(json_res)
    #json = simplejson.dumps(data)

    dict_cur.close()
    data = {'geojson':  mark_safe(json_res)}
    return render(request, 'map.html', data)

urls.py

 url(r'^data/$', GeoJSONLayerView.as_view(model=PFT), name='data'),

models.py

class PFT(models.Model):
    species=models.CharField(max_length=50)
    # GeoDjango-specific: a geometry field (MultiPolygonField)
    geom = models.PointField(srid=4326)
    objects=models.GeoManager()
    # Returns the string representation of the model.
    def __str__(self):              # __unicode__ on Python 2
        return '%s %s %s' % (self.species, self.geom.x, self.geom.y)

How to solve this?


texts: Arcgis IPolygon(System.__comObject) and ILine2 (Esri.Arcgis.LineClass) intersection

texts:
Currently I am developing a Server Object Extension. On my project I am creating a new line class with
LineClass myLine = new LineClass();
and getting a polygon from client as JsonObject than converting this polygon to IPolygon with

JsonObject jsonFieldPolygon;
if (!operationInput.TryGetJsonObject("field", out jsonFieldPolygon))
     throw new ArgumentNullException("field");
IPolygon4 fieldGeom = (IPolygon4)Conversion.ToGeometry(jsonFieldPolygon,esriGeometryType.esriGeometryPolygon);

Polygon I have converted from jsonobject is becoming a COM object, but line I have created is a class in my own application space.

Thus I can't do operations with these two geometries, I guess both must be COM object for operations like intersect and etc.

I am getting an access violation when I try to intersect these two geometries.

My code is below:

JsonObject jsonFieldPolygon;
if (!operationInput.TryGetJsonObject("field", out jsonFieldPolygon))
    throw new ArgumentNullException("field");
IPolygon4 fieldGeom = (IPolygon4)Conversion.ToGeometry(jsonFieldPolygon, esriGeometryType.esriGeometryPolygon);

JsonObject jsonFieldPath;
if (!operationInput.TryGetJsonObject("fieldPath", out jsonFieldPath))
    throw new ArgumentNullException("fieldPath");
IPolyline fieldPathGeom = (IPolyline)Conversion.ToGeometry(jsonFieldPath, esriGeometryType.esriGeometryPolyline);

ITopologicalOperator myOp = (ITopologicalOperator)fieldGeom;

IGeometry resultGeom = myOp.Intersect(fieldPathGeom, esriGeometryDimension.esriGeometry1Dimension); // this line works because both geometries are com objects

IPoint startPoint = new PointClass();
IPoint endPoint = new PointClass();
LineClass myLine = new LineClass();

myLine.SpatialReference = fieldGeom.SpatialReference;
endPoint.SetEmpty();
startPoint.SetEmpty();

startPoint.PutCoords(fieldGeom.Envelope.XMin, fieldGeom.Envelope.YMin);
endPoint.PutCoords(fieldGeom.Envelope.XMax, fieldGeom.Envelope.YMax);

myLine.PutCoords(startPoint, endPoint);

IGeometry resultGeom2 = myOp.Intersect(myLine, esriGeometryDimension.esriGeometry1Dimension);// this line throws accessviolation exception because one geometry is com other is a class in my app

I hope there is a way to solve this


texts:
I am converting multiple CSV (Lon, Lat, value) to GeoTIFF. My output is in a different projection (EPSG:404000) using

output = 'gdal_grid -a invdist:power=2.0 -zfield "value" -of GTiff -ot Float64 -l %s %s %s' % (lyr_name, vrt_fn, out_tif)

I would like to convert it to EPSG:4326 using each of the following lines alone with no luck. Can I add something to my original output to set the correct projection?

output = 'gdaltransform -s_srs EPSG:404000 -t_srs EPSG:4326 -l %s %s %s' % (lyr_name, vrt_fn, out_tif)

outputSRS = 'epsg:4326'

output = gdal_translate -projwin_srs epsg:4326


texts:
I'm creating a bubble map using Leaflet JS. Everything is working smoothly (LeafletJS is awesome!) but now I'd like to try flannery compensation on the bubbles.

I can't find the math anywhere though. Best explainations I've been able to find say you have to multiply radius by either 1.4 or 1.8, which just seems a bit useless.

What is the math used for Flannery compensation?


texts: Listing Attribute Table in order of Selected objects

texts:
Currently when I open my attribute table, I get the order or points by which they are added and I would like to get them in order in which I selected them.

For example,
I select points (in order): 6, 24,3,91
and in my attribute table it shows the points in OBJECT_ID value from lowest to highest.
3,6,24,91

However I would want these to come up in the order I select them. Any ideas how to change this. I'm sure its staring me right in the face.

I have tried to change Advances Sorting to be turned off but this still doesn't fix my issue.


texts: Dropping points based on class number in lidR

texts:
I have a TLS point cloud which were classified into ground (class 2) and unclassified (class 1).
lasinfo (190812) report for 'infile.laz'
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            0
  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000
  version major.minor:        1.2
  system identifier:          'LAStools (c) by rapidlasso GmbH'
  generating software:        'las2las (version 190623)'
  file creation day/year:     50/2019
  header size:                227
  offset to point data:       337
  number var. length records: 1
  point data format:          3
  point data record length:   34
  number of point records:    52463936
  number of points by return: 52463936 0 0 0 0
  scale factor x y z:         0.00025 0.00025 0.00025
  offset x y z:               695760.719696044921875 6106543.521759033203125 0
  min x y z:                  695760.73345 6106543.76051 -0.03700
  max x y z:                  695938.96420 6106715.68576 14.72325
variable length header record 1 of 1:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34735
  length after header  56
  description          'GeoKeyDirectoryTag (mandatory)'
    GeoKeyDirectoryTag version 1.1.0 number of keys 6
      key 1025 tiff_tag_location 0 count 1 value_offset 2 - GTRasterTypeGeoKey: RasterPixelIsPoint
      key 2052 tiff_tag_location 0 count 1 value_offset 9001 - GeogLinearUnitsGeoKey: Linear_Meter
      key 2054 tiff_tag_location 0 count 1 value_offset 9102 - GeogAngularUnitsGeoKey: Angular_Degree
      key 2060 tiff_tag_location 0 count 1 value_offset 9102 - GeogAzimuthUnitsGeoKey: Angular_Degree
      key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
      key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter
LASzip compression (version 3.4r1 c2 50000): POINT10 2 GPSTIME11 2 RGB12 2
reporting minimum and maximum for all LAS point record entries ...
  X                  55     712978
  Y                 955     688656
  Z                -148      58893
  intensity           0      65432
  return_number       1          1
  number_of_returns   1          1
  edge_of_flight_line 0          0
  scan_direction_flag 0          0
  classification      1          2
  scan_angle_rank     0          0
  user_data           0          0
  point_source_ID    64         64
  gps_time 194.000000 21782810.000000
WARNING: range violates GPS week time specified by global encoding bit 0
  Color R 0 61696
        G 0 59904
        B 0 64256
number of first returns:        52463936
number of intermediate returns: 0
number of last returns:         52463936
number of single returns:       52463936
overview over number of returns of given pulse: 52463936 0 0 0 0 0 0
histogram of classification of points:
        34352771  unclassified (1)
        18111165  ground (2)

I want to drop ground class 2 and keep class 1 in analysis:
> las_sf = readLAS("infile.laz", select = "xyzc", filter = "-keep_first -drop_z_below 0 ")

However, it doesn't seem to work:
> npoints(las_sf)
[1] 52463535

I am using lidR_3.0.1
How can I work around this?


texts: Is there a dynamic text method in ArcGIS Pro to mandate the display of trailing zeros in map coordinates?

texts:
I am interested in showing map center coordinates on a layout in ArcGIS Pro using dynamic text:

<dyn type="mapFrame" name="MapName" property="center.y" units="dd" decimalPlaces="4" showDirections = "False"/>, <dyn type="mapFrame" name="MapName" property="center.x" units="dd" decimalPlaces="4" showDirections = "False"/>

This gives me, for example: 36.8101°, -76.612°

This seems to imply that one value is less precise than the other, which shouldn't be the case. I'd like to show the longitude value as "-75.6120" with the trailing 0 (ideally losing both the degree symbols as well, if possible).

I have played around quite a bit with the recommendations provided here, under the Decimal Degree Parts section, without success:
https://pro.arcgis.com/en/pro-app/help/layouts/add-and-modify-dynamic-text.htm

I am likely just missing some combination of tags, but cannot figure it out.


texts: Query SDE to display records created for 'X' date/time

texts:
I am looking to query our SDE to display feature classes modified beginning from January 1, 2013 to today.

I am using Microsoft SQL Server Management Studio, and as it stands searching through dbo tables but can't find any information that would display feature classes modified.

The purpose is to build a business case by identifying a volume of edits/modifications to our data over time.


texts: Need a spatial database that can be used by QGIS and R

texts:
I'm quite new to QGIS and I'm searching  for a simple way to create a Database that allows me to use the tables in QGIS and R. I have a lot of data stored in Excel spreadsheets wich will be used for the statistical analysis in R.

My tasks in more detail:
My first task will be visiualizing the data by joining the tables to a shapefile.
Second I want to make some spatial transformations within QGIS with shapefiles and use the newly created data in my statistical analyse. Therefore I will have to join this data to my database.

Last but not least I wan't to assure that the data is stored in a good way for future analises in QGIS and R since the database will be used in other contexts.

Some tolled me that MS Access is a good way to do so. Can anyone of you give me some advice? Maybe on literature and stuff like that too.


texts: When trying to connect to remote PostGIS database with QGis, getting "ERROR: relation 'geometry_columns' does not exist

texts:
Getting the error message/popup stating:

Connection failed - Check settings and try again.

Extended error information:
ERROR: relation 'geometry_columns' does not exist.

LINE|1: ...attnum as attnum FROM pg_class c, pg_namespace n, geometry_c... ^

when trying to connect to a PostgreSQL/PostGIS database remotely on QGis (Mac). The pg_hba.conf and postgresql.conf files on the server don't appear to be the issue here. Perhaps there's a way I can grant all privileges for these tables mentioned in the error, but there aren't any tables called 'geometry_columns', so I'm stumped so far.

Version of postgresql on remote server is:

SELECT version();
#=> PostgreSQL 9.1.24 on x86_64-unknown-linux-gnu, compiled by gcc (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3, 64-bit


texts:
I have a single table of data (crimes committed) and a region that I want to compare.

I use the Percent Change function of the Crime Analysis Ribbon, I have to feed it two tables and a polygon.

Previously I have exported the table to excel, split it into two tables and re-imported it for the two dates that I want to compare (i.e Q1 vs Q2) however this seems very inefficient and wasteful and I'm sure there's a way in ArcGIS to be able to either do this without splitting the table, or split the table automatically.

The table is super simplistic its just {Lat, Long, Date, CrimeType}


texts: Use pyproj new projection initialization method when setting a CRS in GeoPandas

texts:
Issue

Using GeoPandas:

this:

(...)
t_srs = 4326
gdf.set_geometry('geom', crs={'init': u'epsg:'+str(t_srs)}, inplace=True)

or this:

(...)
gpd.read_postgis(sql_string, conn, geom_col='geom', crs={'init': u'epsg:'+str(t_srs)})

raise a warning today:

/usr/local/lib/python3.6/dist-packages/pyproj/crs.py:77:
FutureWarning:
    '+init=<authority>:<code>' syntax is deprecated.
    '<authority>:<code>'  is the preferred initialization method.
  return _prepare_from_string(" ".join(pjargs))

and if I remove the 'init' key of the crs dictionary, I got this error;

CRSError: Invalid CRS input: {'epsg:4326'}

And if, by chance, I try with {'epsg':'4326'} or {'epsg':4326} ; these errors are respectively raised:

CRSError: Invalid projection: {"epsg": "4326"}: (Internal Proj Error: proj_create: Missing "type" key)

and:

CRSError: Invalid projection: {"epsg": 4326}: (Internal Proj Error: proj_create: Missing "type" key)

Question

Is it already possible to write some GeoPandas code that corresponds to the new* pyproj projection initialization method using  <authority>:<code> instead of  +init=<authority>:<code> ?

https://pyproj4.github.io/pyproj/stable/gotchas.html


texts: Installing MOLUSCE plugin on linux? (Error reading metadata: general)

texts:
I am trying to install MOLUSCE plugin on ubuntu but every time i try i get an error message saying:

The plugin is broken. Python said:
Error Reading metadata: general

There is a screenshot so, does anyone have a solution?
I have used this plugin on Windows earlier but now i am working with linux so i would like to now if there is some solution for using this plugin on ubuntu?


texts: How to sort multi-column lists by first or second column in Google Earth Engine

texts:
I have a list with two-segment.
var list = ee.List([
  [101, 102],
  [102, 103],
  [103, 104],
  [106, 100],
  [104, 105],
  [105, 106],
  ])

How can sort this list by first or second column?


texts: Select intersecting features within Value Relation Widget

texts:
Is it possible to filter records in the Value Relation Widget to include only those that intersect the current feature?

We are trying to set up a layer for user input where we need to populate Start/End cross-street intersections. I can use the Value-Relation widget to bring up a list of all the Street-Names from the master-roads layer, but I'd like to reduce it to jus those roads that intersect the feature being edited. As per this example, sometimes the road changes name at either side of the intersection, so for consistency it would be best if manually selected by the data-entry person rather than automatically populating and just ending up with the first match (or whatever a computer decided!).

I'm not sure if the filter expression can be used in this way, but it'd be great if it can. I've tried a few expressions but so far no success.


texts: GeoWebCache + OpenLayers - caching blank tiles?

texts:
Are browsers supposed to cache blank tiles? I have a layer cached with GeoWebCache. When I request this layer via OpenLayers I can see that tiles containing data are taken from cache, the blank tiles are always requested from GeoServer (and then taken from cache?). Here is the screen - all requests on first screen returned empty tiles, still it takes lot of ms:

Here are the tiles containing actual data:

Is it standard behaviour of OpenLayers or I did something wrong? (I am using OpenLayers 2.13.1).

I can't find any header that may be responsible for caching by browser. Is it possible that OpenLayers is doing the caching based on the response?

Empty tiles response:

Tile with data response (cached)

So now I have 3 new questions:

What is the reason? Is it because GeoWebCache is not caching empty tiles? Or OpenLayers is not letting browser to cache them?
Is it more optimal not to cache empty tiles?
How to change this behaviour?


texts: I can't figure out why the geom field is so long, I can't make heads or tails of the number either.


texts: Loading Open Data Cube data to match a raster with geographic coordinates?

texts:
I would like to load a raster from a GeoTIFF file, and then load a time series of data from the Open Data Cube into the same CRS and resolution so I can use the raster as a mask. I can use the raster's .geobox to do this nicely when my GeoTIFF is in projected coordinates:
import xarray as xr
import datacube
dc = datacube.Datacube()

# Load raster and remove the redundant "band" dimension
raster = xr.open_rasterio("raster.tif").squeeze("band")

# Load data from datacube into raster geobox
ds = dc.load(
    product="s2b_ard_granule",
    measurements=["nbart_red"],
    time=("2018-01", "2018-02"),
    like=raster.geobox
)

# Use raster data as a mask
ds.where(raster > 500)

However, if my raster is in geographic coordinates (degrees), datacube will load with dimensions named longitude and latitude, not x and y. This means the dimensions no longer match the dimensions in raster, which prevents me from using it as a mask in the final step.
I've tried simply re-naming the dimensions in ds using ds.rename({'longitude': 'x', 'latitude': 'y'}). This fixes the immediate problem and allows me to use the raster data as a mask, but ds  still has a .geobox with dimensions named longitude and latitude after the rename which now don't match the xr.DataArray's dimensions. This causes further issues downstream where I want to use the .geobox for other steps.
Is there:

Any way to load data from datacube to match a raster .geobox, but not switch to the longitude and latitude  dimension names even if the raster has geographic coords?
If not, is there any way to completely rename the coords so that they are renamed in both the xr.DataArray and the .geobox?


texts: RowChangedEvent doesn't catch every change

texts:
I'm using C# to create ArcGIS Pro add-ins.

Until now I've been able to catch changes made directly on the attribute table. This is modify, create or remove a row. But, when I use a geoprocessing tool like 'Calculate Field", my code won't catch this change.

This is my script now:

using ArcGIS.Core.Data;
using ArcGIS.Core.Events;
using ArcGIS.Desktop.Editing;
using ArcGIS.Desktop.Editing.Attributes;
using ArcGIS.Desktop.Editing.Events;
using ArcGIS.Desktop.Framework.Contracts;
using ArcGIS.Desktop.Framework.Dialogs;
using ArcGIS.Desktop.Framework.Threading.Tasks;
using ArcGIS.Desktop.Mapping;
using System;
using System.Collections.Generic;
using System.Linq;

namespace AddinsArcGISPro
{
    internal class MonitoringButton: Button
    {
        private static SubscriptionToken _rowChangedToken;

        protected override void OnClick()
        {
            var featureLayer = MapView.Active.Map.Layers.First(t => t.Name == "TableName") as FeatureLayer;

            QueuedTask.Run(() =>
            {
                _rowChangedToken = RowChangedEvent.Subscribe(OnRowChangedEvent, featureLayer.GetTable());
            });
        }

        private static void OnRowChangedEvent(RowChangedEventArgs obj)
        {
            // Do something
        }
    }
}

What should I do to catch those changes made by geoprocessing tools?


texts: config parser error raised by TileCache and Mapserver with OpenLayers

texts:
Just a forewarning, I will try to be as extensive with my question as possible. Due to a number of reasons why this might not be working, I have posted all sorts of portions of my debugging below.

Web map located here

Tile cache is being setup as a test run for the "hii_5" layer in OpenLayers. MapServer is the backend.

I am trying to get TileCache to cache images of my MapServer WMS. However, no images are produced. The TileCache is "trying" to load the tiles (as you can see on the live page), however, no images are returned.

When I follow the path of one of the image tiles, I see this error:

An error occurred: File contains parsing errors: C:\ms4w\Apache\cgi-bin\tilecache.cfg
    [line 82]: '  type=WMS\n'
    [line 83]: '  url=http://141.117.104.183/cgi-bin/mapserv.exe\n'
    [line 84]: '  mapfile=C:/ms4w/Apache/htdocs/moot2_dev/hii_railroads.map\n'
    [line 85]: '  extension=png\n'
    [line 86]: '  layers=HII_railroads\n'
  File "C:\ms4w\Apache\cgi-bin\TileCache\Service.py", line 85, in _load
    config.read(files)
  File "C:\Python27\lib\ConfigParser.py", line 297, in read
    self._read(fp, filename)
  File "C:\Python27\lib\ConfigParser.py", line 538, in _read
    raise e

I've run a GetCapabilities request and there are no warnings
but GetMap request returns a blank page
and the template return blank

Here are the relevant mapfile, tilecache.cfg, and openlayers source:

Not sure how much help I will get here, but the questions out there for anyone to take a stab at it!

hii_railroads.map

MAP                                                     #-> BEGIN MAP
    OUTPUTFORMAT                                        #-> BEGIN OUTPUT FORMATTING DEFINITION
        NAME "png"                                      # Name of format
        DRIVER AGG/PNG                                  # Image driver
        MIMETYPE "image/png"                            # Type of image
        IMAGEMODE RGBA                                  # Image mode
        EXTENSION "png"                                 # Image filetype
        TRANSPARENT ON                                  # Transparency
    END                                                 #-> END OUTPUT FORMATTING DEFINITION

    WEB                                                 #-> BEGIN WEB DEFINITION

        METADATA
            "wms_title"           "HII_railroads"
            "wms_onlineresource"  "http://rinnerweb.arts.ryerson.ca/cgi-bin/mapserv.exe?map=/ms4w/Apache/htdocs/moot2_dev/hii_railroads.map&"
            "wms_srs"             "EPSG:3857"
            "ows_enable_request"  "*"
        END

        TEMPLATE "template.html"
        IMAGEPATH "/ms4w/tmp/ms_tmp/"
        IMAGEURL "/ms_tmp/"
    END                                                 #-> END WEB DEFINITION

LAYER                                                   #-> BEGIN LAYER DEFINITION: Railroads <-#

    METADATA
        "wms_title" "HII_railroads"
    END

    NAME    HII_railroads                               # Layer identifier
    DATA    "data/wgs84/railroads_360.tif"              # Name of the data (in this case, string identifying a raster)
    STATUS  ON

tilecache.cfg

[HII_railroads]
  type=WMS
  url=http://141.117.104.183/cgi-bin/mapserv.exe
  mapfile=C:/ms4w/Apache/htdocs/moot2_dev/hii_railroads.map
  extension=png
  layers=HII_railroads

index.html

hii_5 = new OpenLayers.Layer.WMS("Railroads","../../cgi-bin/tilecache.cgi?",{
    layers:"HII_railroads",
    format:"image/png"},{
    isBaseLayer:false,
    transparent:true,
    format:"image/png",
    alpha:true,
    displayInLayerSwitcher:false
    });


texts: Access specific version of featureClass (Table) through SOE

texts:
I am new to ArcObject and SOE.

The task I am trying to accomplish is:

Find the specific version of a table
Copy the table (to a new table) that the non-soe user can access that version of the table.

Any suggestions with which direction I should go?


texts: Automated way of Connecting dangling polylines of two separate feature classes

texts:
I have two datasets that I am trying to form into a network, one representing pedestrian and bike paths, the other representing roads.

The lines do not connect together (as seen in the picture below) and I would like to know if there are any automated ways to connect the end-points of each line segment from the two layers together?


texts: Finding new point, with proportion, between two points Postgis

texts:
I couldn't find a topic related to my problem, thus I am going to ask. I have the following problem:
I need to create a new Point based on three variables: P1 (Geometry), P2 (Geometry) and length (Integer in meters). I know the geometry of two points and want to crate a new point the lies between them given a determined length, like in the picture below:

How can I do it using only Postgis?
I am using Postgis 2.2 with PostgreSQL 9.4.


texts:
I'd like to create custom buffers in QGIS. Each buffer will surround a point. The buffers don't overlap. And each buffer has to have two constants:

There is a 50m radius space around the point which is within the buffer.
The buffer as a whole covers 10 acres.

The shape of the buffer will vary from buffer to buffer, so this detail is immaterial. What's important is, the buffers cover ten acres, including the 50m.

The above illustration is an example of what I'd want the polygons to look like, they can be of any shape in the world - so long as  each one covers, in total, 10 acres.


texts:
I made a leaflet map with a search control. The problem occurs when enabling circleLocation and flyTo: The circle gets drawn much too early and in the center of the current mapbound. When the map reaches the position (via "flyTo") of the point i was searching for, only then the circle switches to his correct position (around the location found). To work corretly the circle should only be drawn after flyTo is finished.

I already tried it with panTo or setView and everything worked fine. But flyTo looks so much better and i would like to use it.

Thats my code:

    showLocation: function(latlng, title) {
    if(this.options.zoom)
        this._map.setView(latlng, this.options.zoom);

        this._map.flyTo([mycoordinate, mycoodinate2], 11.5);

    if(this._markerLoc)
    {
        //clearTimeout(this._markerLoc);
        var that = this;
        this._markerLoc.setLatLng(latlng);
        this._markerLoc.setTitle(title);
        this._markerLoc.show(), setTimeout(function() {
        that._markerLoc.hide();
        }, 1250);

        if(this.options.animateLocation)
            this._markerLoc.animate();

    }

    if(this.options.autoCollapse)
        this.collapse();
    return this;
}});


texts: Using ModelBuilder and Inline substitution in ArcGIS Pro?

texts: I am creating a model and was wondering whether it is possible to use inline substitution with variables that are a table of values or multiple values, shown below.


texts:
It's my understanding that in 2008, the US government published GPS standards for civilian applications (including cellphones) of rms 4m.

Is it correct that this applies to cellphones? and
Has this 4m rms been improved since 2008?

http://www.gps.gov/technical/ps/2008-SPS-performance-standard.pdf

What is the maximum Theoretical accuracy of GPS?


texts:
We have a maritime simulator in our office we occasionally use for engineering studies...I was trying to figure out how to get a shape file of the path taken by the ship being driven into QGIS.  Data is being sent to a PC with a charting software called Rosepoint on it in NMEA format over a TCP stream via a network cable.  According to Rosepoint it is doing this via NMEA over TCP.
Initially hooking up QGIS I didn't think it was seeing anything, I tried Global Mapper which really saw nothing as it wouldn't even connect.  It did however occur to me that with QGIS my lat/long were both 0 so it must at least be connecting with autoconnect.  If I set the com port to something else it fails to connect all together.
So I'm assuming maybe this is working and it's just not able to understand the sentences coming through...is there a way to define the format on my own, or something I'm not doing correctly??


texts: qgis2threejs viewer showing poor resolution

texts: When I open qgis2threejs in preview mode the resolution of the preview is poor. It used to be very good and would render as you moved around eth preview screen.


texts: Displaying dynamically changing data layer in moving view using OL-Cesium

texts:
I'm developing an application that contains a tilted map based on OL-Cesium. The "camera" height is very low (zoom level is the highest), and the view is "moving" continuously, like a car navigation in 3D view. The moving speed is between 0 and 70 km/h. The map has some tileWMS layers published by GeoServer. There are no problems with them. However, there is a polygon layer with dynamically changing content (inserts in about every second). The source of this layer is a PostGIS table.
Now, I'm trying some methods to put this data onto the map. I think that the tileWMS is not an appropriate choice for dynamically changing data. Now, I tried the ImageWMS. It performs well in 2D view, but in 3D it is very slow (and the polygons look very bad due to the high scale).
The second option is the WFS with GeoJSON format and BBOX loading strategy. In this case, the problem is the very high CPU and the RAM using is growing continuously. So, now I'm trying to find a better solution:

Call a function with a defined frequency (50-100 ms) which selects the polygons from the PostGIS table that are within X meters like this:
SELECT ST_AsText(ST_Collect(geom)) FROM table t WHERE
ST_Dwithin(t.geom,ST_SetSRID(ST_Point(x_coord,y_coord),4326)),0.003);

Store the results in WKT format in a string.

Pass the WKT string to the OpenLayers.

Refresh the data source of the vector layer like this:
function refreshVectorSource(polygonWKT) {
  vectorSource.clear();
  var format = new ol.format.GeoJSON({
    featureProjection: "EPSG:4326"
  });
  vectorSource.addFeatures(format.readFeatures(polygonWKT))
};

What do you think? Do you have any suggestions? Or any other idea?


texts: How to retrieve current extents of the location?

texts: How google is getting its maps zoom to current extents where we have located eg: If i am get hyderabad kphp on google maps its directly redirects to hyderabad city extents while opening a map , if it is chennai the if i open google maps in chennai the city extents automatically changes do it goes ? Is it based on IP address location or what?


texts: Is MULTICURVE geometry in WKT valid according to the standard and open source community?

texts:
We are using Oracle function (SDO_UTIL.TO_WKTGEOMETRY) for exporting geometry to WKT. Some geometries have MULTICURVE in WKT, like this one:

MULTICURVE (COMPOUNDCURVE (LINESTRING (3451418.006 5481808.951, 3451417.787 5481809.927, 3451409.995 5481806.744), LINESTRING (3451409.995 5481806.744, 3451418.006 5481808.951)), LINESTRING (3451418.006 5481808.951, 3451417.787 5481809.927, 3451409.995 5481806.744), LINESTRING (3451409.995 5481806.744, 3451418.006 5481808.951))

Then we store the geometries in PostgeSQL with PostGIS extension and geometry is valid from PostGIS point of view. But then we are using WKB parser in our app which does not support parsing for geometries with MULTICURVE instance. In the implementation of simple features we found on the page 24 that MULTICURVE is only abstract class and instances shouldn't be created.

But Oracle export geometry with MULTICURVE. PostGIS reads the MULTICURVE and QGIS also. We found some examples on the internet in the opensource libraries which are able to parse MULTICURVE geometries and are inluded in the examples of WKT.

Now we don't know if we have a bad standard or open source community use MULTICURVE despite of standard or if there is something else which we don't know.


texts: Can't insert AI CS6 files as objects into ArcMap 10.1

texts:
In the past I have inserted Adobe Illustrator files into ArcMap as objects. My old version of AI was 8.0 and ArcMap was 9.2. Now I have arcmap 10.1 and AI CS6. When the process doesn't shut down arcmap, it gives the message.."failed to create object. Make sure object is entered in system registry.
???????


texts: Working on ARCGIS for Offline mapping ,Can a tpk file have multiple map layers inside it?


texts: AutoCAD Map/Civil 3D Feature Data Object (FDO) layer displays but is no longer editable

texts:
Just over a week ago, one of our main GIS layers became uneditable. Our GIS layers are connected to AutoCAD Civil 3D (Fig.3) via an SQL Server Spatial Connection in the Task Pane. We were able to edit all connected layers for years, until last week when Cadastre became read only within AutoCAD Civil 3D. All other Feature Data Object (FDO) layers remained editable (Fig.1), while Cadastre is only visible and selectable (Fig.2).

What we know up to this point:

Cadastre is selectable but cannot be edited in AutoCAD Civil 3D
All other layers can be selected and edited in AutoCAD Civil 3D
DB Admin reported no locks or unusual activity on the Cadastre table
All layers, including Cadastre, are editable in another application - QGIS
Cadastre table was reconnected several times in AutoCAD Civil 3D
This affects all users using AutoCAD Civil 3D on separate workstations
All users have the same edit permissions to the Cadastre table in the database

Can you help us resolve this issue?

Fig.1 - Other layer objects present edit grips when selecting them.

Fig.2 - Cadastre can be selected but no edit grips presented. No edit commands or function are available for this layer.

Fig.3 - AutoCAD Civil 3D version.


texts: Definition Query for records created in the past 30 days in and SQL SDE?

texts:
I am stumped as to how to get this to work in ArcMap...

I have records being created ever day, and would like to be able to show records on a map that have been created within the last 2 days based on a date field in our SDE feature class. How can I write a definition query in SQL to do this so that it automatically changes to the current date and then grabs records from the past 2 days? I've looked on the web, but have had no luck finding anything that matches what I need.

The values in our SQL 2012 database for the created_date field are formatted as such:

2014-11-04 18:44:09.0000000


texts: How to create a QgsFeature with default attributes?

texts:
I have a polygon layer loaded from POSTGIS database. This layer has a default value for the first column. I can create a new feature set that attribute to some predefined value with

    newFeature = QgsFeature()
    newFeature.setFields(activeLayer.pendingFields())
    newFeature.setAttribute(0, "someValue")

But I want to set it to the default value that comes from the layer. When I try to add a feature with QGIS 2.10 it automatically sets the default value for that column. How can I achieve that with PyQGIS?


texts: gdal RasterizeLayer gives different extent to arcpy PolygonToRaster

texts:
I have an existing ArcPy script which uses PolygonToRaster_conversion to convert a polygon shapefile to raster (tiff). I am now exploring GDAL/OGR and have used the cookbook approach to carry out the same task but find two key differences in the output rasters:

The extent of the GDAL output raster is different from the arcpy output. Xmin is the  same in each. However, the number of rows and cols are one less each in the gdal rasters, and the Y values differ by an odd amount (which is not a multiple of grid pixel size). I have tried the rasterize function with and without the ALL_TOUCHED option to no avail. I suspectthe difference is in calculating the extent, and applying the GeoTransform. Any clues?
The NoData value set using gdal is neither correctly set or recognised in ArcGIS or QGIS. Instead of a NoData value, the pixels have a value of zero. Could this be to do with the output raster datatype?

My code using gdal:

    NoDataVal = -9999

    # Open the data source and read in the extent
    shpDS = ogr.Open(inPolygonShp)
    shpLayer = shpDS.GetLayer()

    # Create the destination data source
    xMin, xMax, yMin, yMax = shpLayer.GetExtent()
    xRes = int((xMax - xMin) / inGridSize)
    yRes = int((yMax - yMin) / inGridSize)
    rasterDS = gdal.GetDriverByName('GTiff').Create(outputRaster, xRes, yRes, 1, gdal.GDT_Byte)

    # Define spatial reference
    rasterDS.SetProjection(shpLayer.GetSpatialRef().ExportToWkt())
    rasterDS.SetGeoTransform((xMin, inGridSize, 0, yMax, 0, -inGridSize))

    rBand = rasterDS.GetRasterBand(1)
    rBand.SetNoDataValue(NoDataVal)
    rBand.Fill(NoDataVal)

    # Rasterize
    err = gdal.RasterizeLayer(rasterDS, [1], shpLayer, burn_values=[1], options = ["ALL_TOUCHED=TRUE"])     # http://www.gdal.org/gdal_rasterize.html

and using arcpy:

    idFld = arcpy.Describe(inPolygonShp).OIDFieldName
    addMsgAndPrint("OID Field used: "+ idFld)
    arcpy.PolygonToRaster_conversion(inPolygonShp, idFld, outputRaster, "CELL_CENTER", "NONE", inGridSize)


texts: Using ArcGIS API for Python and Windows Authentication?

texts:
I've been trying to follow the ESRI recommended workflow  to log on to Enterprise (using authentication), but it just doesn't want to work. I'm able to get the account credentials registered in Windows Credential Manager, but if I try to run the program and access them via the API for Python, I keep receiving this error:
AttributeError                            Traceback (most recent call last)

<ipython-input-1-9fd4325b765b> in <module>

     11

     12

---> 13 gis =GIS(portal_url,profile=portal_profile, verify_cert=portal_cert_check)

     14 print('Good')

~\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone1\lib\site-packages\arcgis\gis\__init__.py in __init__(self, url, username, password, key_file, cert_file, verify_cert, set_active, client_id, profile, **kwargs)

    337                                    "argument when connecting to the GIS.")

    338             else:

--> 339                 raise e

    340         try:

    341             if url.lower().find("arcgis.com") > -1 and \

~\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone1\lib\site-packages\arcgis\gis\__init__.py in __init__(self, url, username, password, key_file, cert_file, verify_cert, set_active, client_id, profile, **kwargs)

    321                                            expiration=self._expiration,

    322                                            referer=self._referer,

--> 323                                            custom_auth=custom_auth)

    324             if self._is_hosted_nb_home:

    325                 # For GIS("home") objects, force no referer passed in

~\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone1\lib\site-packages\arcgis\gis\_impl\_portalpy.py in __init__(self, url, username, password, key_file, cert_file, expiration, referer, proxy_host, proxy_port, connection, workdir, tokenurl, verify_cert, client_id, custom_auth, token)

    165                                       client_id=client_id,

    166                                       custom_auth=custom_auth,

--> 167                                       token=token)

    168         #self.get_version(True)

    169         self.get_properties(True)

~\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone1\lib\site-packages\arcgis\gis\_impl\_con\_connection.py in __init__(self, baseurl, username, password, **kwargs)

    156             self._auth = "PKI"

    157

--> 158         self._create_session()

    159

    160         #  Product Info

~\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone1\lib\site-packages\arcgis\gis\_impl\_con\_connection.py in _create_session(self)

    260             self._session.auth = GuessAuth(username=self._username,

    261                                            password=self._password)

--> 262         elif self._auth.lower() == "pro":

    263             self._session.auth = GuessAuth(None, None)

    264         else:

AttributeError: 'NoneType' object has no attribute 'lower'


texts: Selecting highest value within group using QGIS?

texts:
I have a .shp of points. This points are organized in groups with same ID (group 1, 2 ,3)

Each group has 10 points.

Each point has a Z column with Z values.

I want to select the point with the highest Z value for each group.

So just one point should be selected for each group.

I was trying "select by expression" but I don't know what expression to use.

Also, I need a way to automatize this. This process is just a small part of my project that needs to be repeated a bunch of times, so doing this manually all the time is impossible. A way to use do this with graphical modeler would be great.


texts: How to set an x-offset for my whole linestring/layer?

texts:
Supposing I have a layer imported from GeoServer succesfully, and it looks like this:
(some will remember me from my previous question. And, yes, that is my city Volos, Greece. There is a map and a blue MultiLinestring)

Supposing this layer is offset by a few pixels by default to the left, which is the way to put it back where I want it.

How can I offset the whole LineString a few pixels to the right?

I tried this way:

var offsetStyle = new OpenLayers.Style({
    graphicXOffset: 50000,
    graphicYOffset: 50000
});

wms.style = offsetStyle;

But it doesn't work! The LineString will not move to the left or the right!

Can I do anything with QGIS to move my layer to the left or the right?


texts: I'd like to print out 6 images in a single page. How can I do this using map composition tools in ENVI (6 images are year-to-year dynamic images  with the same size and the same geo-referencing)?


texts:
When I harvest my GeoServer WFS using GeoNetwork to create iso19139 records, the URL of the first entry within MD_DigitalTransferOptions is always:

http://localhost:8000/geoserver/wfs?request=GetFeature

This isn't particularly useful as that URL doesn't reference the feature and therefore can't be used as a service for that resource. The remaining entries in MD_DigitalTransferOptions correspond to different WFS outputFormats and do specify the typename within the URL element e.g:

http://localhost:8000/geoserver/wfs?request=GetFeature&request=GetFeature&service=WFS&typename=sf:archsites&outputFormat=KML&version=1.0.0

Is there a way to prevent the harvester creating the first entry e.g. either within GeoNetwork or possibly GeoServer WFS settings?