texts: Worldwide High-Resolution Water Body Shapefiles with Consistent Scale

texts:
Where from your website can I access worldwide vector water body shapefiles that have better resolution than the generic SRTM 90m data?

Do you have coastlines and inland features (lakes, estuaries, reservoirs, and rivers) which are better than SRTM and that the scale is consistent across the entire coastline?

Also, is your data surveyed, generated from an open source mapping team, or is it gathered from other sources such as the USGS, Geoscience Australia, etc?

If your data is gathered and collected, do you have published metadata associated per set?


texts: QGIS Spatial Join Includes Records From Features Just Outside Polygon

texts:
I am using QGIS to do spatial joins between UK parliamentary constituency boundaries and subsidy amounts associated with areas defined by boundaries in a separate layer. The idea is to find the subsidy amounts for the parliamentary constituencies. What I am finding is that the join is including areas that appear on screen to be outside the constituency boundary, but with part of the boundary very close to the constituency boundary. My subsidy data is for England only, but the join is assigning subsidies in English constituencies to the neighbouring Welsh constituencies.

Is this a common problem? What is the likely cause? How can it be dealt with? Any help very welcome.


texts: Is it possible to (trivially, via a tool in python) to convert a polygon to a linestring? if so, how? If not, why?


texts: Keeping same color when using "gdalogr:polygonize" function to polygonize segmented image

texts:
This is one example of images that I have:

Now I want to polygonize it to get shapefile of different polygons. Every time I run this algorithm I get polygons of different colors but I want to keep colors of polygons same as on image parts that is made from.  Is it  possible, does anyone have code for this?


texts: How to set a maximum amount of features to be drawn in every zoom level

texts:
I want to apply a form of generalization/specialization to my OpenLayers map. Is there a way to set the amount of WFS features to be rendered based on the zoom level?

In GeoServer you can set the maximum amount of features but not based on a zoom level. Would something like this be possible?


texts: How to get a drawn feature id in OpenLayers

texts:
Using OpenLayers 4. I can get the list of drawn features by calling getFeatures(), but calling getId() on each feature returns "undefined".
Here is a screenshot from console with the list of drawn features.

As far as I see "c" is the "id" which is returned by getId() and  is "undefined" for all drawn features by default. So do drawn features have some automatically generated unique ids? The "xp" value seems to be the one. how can I get it then?


texts: Subtracting rasters from each other in batch using ArcGIS for Desktop?

texts:
I am trying to calculate the maximum change for each cell in a large amount of rasters. These rasters are number X002,X003 till X150 and I am trying to calculate the change between two consecutive rasters and eventually checking the maximum change per cell by comparing all the changes in the raster values.

I have tried doing this by using a raster iterator to load the rasters into the model, this is the point I get stuck at. The problem is that I can't get to load both raster X002 and X003 at the same time and thus can't calculate the change.


texts: Can a shapefile have polygons that overlap each other?

texts:
As the title says: a simple question!

Can a valid shapefile have polygons that overlap each other? Or would they need to be in separate shapefiles?

I have Googled for the answer, but with no success...

UPDATE: Found the shapefile technical description, which says "Shapefiles handle single features that overlap or that are noncontiguous". Though I find that ambiguous as to whether that means a single self-overlapping feature, or single features that overlap each other, so my question stands.


texts: How can I do multiple snapping of vector layers in QGIS?

texts:
I have a network of a transport system that I created from 1970-2013. However, the accuracy is letting me down as I used different snapping settings throughout the process. Is there a plugin or a way that I can correct each line (and align each method of transport which uses the same stretch of road) without using the node tool?

I have attempted to use the traceDIGITISE plugin but it didnt work for me.


texts:
I just started exploring the CMV about a week ago as an option to replace our Flex Viewer based configurable web mapping applications. So far it seems to be meeting all our needs however I have run into a snag regarding layer order.

Our web maps consume a variety of operational layer types but when I add a KML layer, it gets placed on top of all other layers. All layers added from our GIS servers get placed in the same order as listed in the code of viewer.js but the KML layers are not honored by this order. Is there any way to manage the layer order of KML layers directly in the CMV?

To clarify, when GIS server layers are mixed with kml layers, the kml layers are grouped separately from the GIS server group of layers regardless of operational layer order in viewer.js. Although, within both of these "groups," layer order as per the code of viewer.js is honored. It's in the rendering of layers on the map that layer order is not honored within the "group" of kml layers. The rendering is completely random, see below for three separate renderings of a map with all the same settings.


texts: Triggering event from marker by pressing enter when marker has focus using Leaflet js

texts:
I need to create a leaflet.js map within certain accessibility guidelines.

1) I would like to be able to trigger an event (not just open a popup) by pressing 'enter' when a marker has focus.

2) I would like to be able to trigger an event on focus of a marker

Is any of this possible in leaflet... I have searched a lot but can't find anything on this.

Some basic pseudo code below for illustration purposes -

 addMapMarkers () => {

   this.mapMarkers.map((item) => {

      const mapMarker = L.marker(this.rc.unproject([item.coordinates.latitude, item.coordinates.longitude])).addTo(this.layerGroup)

      L.DomEvent.addListener(mapMarker, 'click', () => {
        // 1) This works fine obviously - BUT I would also like it to trigger when pressing enter when the marker has focus
        this.map.panTo(mapMarker.getLatLng())
      })

      L.DomEvent.addListener(mapMarker, 'focus', () => {
        // 2) For example pan to center when marker has focus
        this.map.panTo(mapMarker.getLatLng())
      })

    })

  this.map.addLayer(this.layerGroup)
}

The map markers in Leaflet are given a tabindex="0" so as one tabs though a page, on reaching the map one can navigate with the keyboard to zoom etc. On continuing to tab the first marker will come into 'focus' then the next and so on. This suggested to me that 'click' behavior should be available by pressing enter. Google chrome certainly returns a marker as having focus.

document.activeElement

<img src="https://unpkg.com/leaflet@1.5.1/dist/images/marker-icon.png" class="leaflet-marker-icon leaflet-zoom-animated leaflet-interactive" alt="" tabindex="0" style="margin-left: -12px; margin-top: -41px; width: 25px; height: 41px; transform: translate3d(94px, 721px, 0px); z-index: 721; outline: none;">

In response to TomazicM comment....
It appears that clicking enter when a marker has focus will open a popup if present. It does not however trigger an onclick event...

For example:

      L.DomEvent.addListener(mapMarker, 'click', () => {
        // 1) This works fine obviously - BUT I would also like it to trigger when pressing enter when the marker has focus
        this.map.panTo(mapMarker.getLatLng())
      })

Am I missing something here..?

RE: The pull request that was not me.


texts: I am using ArcGIS viewer for Flex 2.4. I am trying to get the overview map to print as part of the map using the print widget. In my customized print widget, the user has the facility to opt for overview map. For this, I intend to open the overview map widget from Print widget. Since it is a UI widget and does not seem to entertain widget id, how do I go about it?


texts: Layer Preview OpenLayers not working: Scale = 1:NaN

texts:
I am taking the tutorial from this site:
http://workshops.boundlessgeo.com/tutorial-autocomplete/

After I published taxlots layer and went to Layer Preview - Openlayers, the map was blank and on the bottom of the box it showed 'Scale = 1 : NaN'. If I chose Preview as JPEG the map showed up fine. It seemed there's something to do with map bounding or reference system? I have doubled checked everything but didn't see anything wrong.

I followed the same procedure and published another layer in shapefile format and in a different coordinate system. It worked fine.

BTW I am using OpenGeo Suite 4.7 and JRE 7 on Windows 7.


texts:
I've been struggling with this problem for awhile. Here is a curvature raster to illustrate the problem.

I'm using 10m DEM rasters from the NRCS Data Gateway for a large portion of California. The problem is some strange tiling that appears after merging, but isn't present in the original data (before merging). The grids create low points, so flow accumulation etc is drawn into them.

I've tried mosaic to new raster, create raster catalog and filter in ArcGIS 10.2, merge and build vrt in Qgis2.4 but the problem persists throughout.

Something to note is that there are apparently two sets of grids, one with regular equant shapes, and then one set where tiles overlap (on both edges of the tile). I discovered this with "build raster catalog" which creates a feature displaying the edges of the files in the catalog.

The lines are flat lows, which creates a problem trying to do any sort of slope stability or hydrologic analysis...which is the whole point of this exercise.

use BILINEAR interpolation (in environment settings) for DEM, these grids are caused by the default nearest neighbor interpolation.


texts: Creating GRASS scripts to run in C# on Windows?

texts:
Most of the tutorials that I have seen are under Python and Linux.

Is it Possible to create GRASS scripts in C# or atleast in Windows?

I can create scripts using a different language as a compromise but at least it must run in C# because I am currently using C# for my GUI of the simulation system that I am making.


texts:
I am trying to do some Line of Sight analysis on some points using 3D Analyst in ArcGIS Desktop. I have a set of points that I've made two copies of a feature class of, one to use as the observers and one to use as the targets (I want to look at line of sight between all the points looking at each other). I made sight lines between them and ran the LoS. The process worked fine. Here's the problem I'm having.
I have elevation data underneath as the raster layer, and am trying to use an OffsetA and OffsetB of 10m to replicate antennas at each point. When I build my sight lines and look at the height of those lines at the beginning and end, it reflects the 10m offsets. But when I run my LoS and look at the Z at the ends of those lines, it shows the height as the elevation data height and doesn't add in the offset.
I've tried using both projected and not-projected. I also tried making my sight lines 2D and 3D. Nothing makes the LoS lines raise up those 10m on the ends. I'm wondering now if those lines are supposed to show at a height with the offset like the sight lines or if they just show from the ground level?
Anyone have any ideas?


texts: Quick OSM request returns a layer which does not match OSM in QGIS

texts:
I made this request from quick osm
[timeout:2500]
[maxsize:2000000000];
area[admin_level='6'][ref='90']->.area;
rel(pivot.area)->.rel;
(
node[highway](area.area);
node[highway](r.rel);
way[highway](area.area);
way[highway](r.rel);
rel[highway](area.area);
rel[highway](r.rel);
);
(._;>;);
out;

The request works fine but when I want to put the layer of Open Street Map, both layers don't fit together. I tried to change the SCR and when the layer from Quick OSMis in ESPG 4326 and Open street map in 2154 they match.
How can I fix that? My QGIS version: 3.10.14


texts: how to visualize csv data every 10 rows in QGIS

texts:
I have imported a csv file in QGIS. The file contains 1 column with 43000 rows of data. I want to apply a filter to visualise data only every 10 rows in order to make more clear the graphics.

Please can you suggest me which type of filter I can use ?


texts: Port Number 8080 required for GeoServer?

texts: Is it necessary to use the port number such as 8080 for GeoServer i.e. will it not work with a different port number?


texts: How do I fix my SLD tile heat map when it comes up as dots?

texts:
I am currently working on an SLD heat map on GeoServer that uses tiles.There is a point in each tile with the data. From far away, the heat map looks normal.

But up close it just becomes little dots within each tile

How do I fix this?


texts: How to Map Latitude/Longitude to FAA Chart?

texts:
How can I map latitude/longitude to an FAA digital chart? Attached is a sample GEOTIFF. More specifically, how can I get a reference points from the digital data?

Once I can get reference points, I have no problem translating Pixel-Lat/Long. In my testing, I have been manually identifying the reference points to the the translation. This is a rather tedious process.

It looks like I just need a single lat/long-pixel-x/y mapping to derive them all from this data.

The GEOTIFF spec suggests the Model Tie Point is the key but the example in the spec does not match this data. -822739.829757951 -64249.12941973 does not give units not does it give offset from what?

In this chart the false origin is located at the Pixel coordinates (6521 3865) from the upper left corner.

Thanks

ExifTool Version Number         : 9.53
File Name                       : CF-19 WAC 43 South.tif
Directory                       : .
File Size                       : 47 MB
File Modification Date/Time     : 2013:03:15 10:45:54-04:00
File Access Date/Time           : 2014:02:28 13:45:59-05:00
File Inode Change Date/Time     : 2014:02:28 13:06:22-05:00
File Permissions                : rwxr-xr-x
File Type                       : TIFF
MIME Type                       : image/tiff
Exif Byte Order                 : Little-endian (Intel, II)
Image Width                     : 11888
Image Height                    : 4135
Bits Per Sample                 : 8
Compression                     : Uncompressed
Photometric Interpretation      : RGB Palette
Strip Offsets                   : (Binary data 286 bytes, use -b option to extract)
Samples Per Pixel               : 1
Rows Per Strip                  : 128
Strip Byte Counts               : (Binary data 262 bytes, use -b option to extract)
Planar Configuration            : Chunky
Color Map                       : (Binary data 1536 bytes, use -b option to extract)
Sample Format                   : Unsigned
Pixel Scale                     : 127.09668 127.09668 0
Model Tie Point                 : 0 0 0 -822739.829757951 -64249.12941973 0
GDAL Metadata                   : <GDALMetadata>.  <Item name="PyramidResamplingType" domain="ESRI">NEAREST</Item>.  <Item name="RepresentationType" sample="0">THEMATIC</Item>.  <Item name="STATISTICS_MINIMUM" sample="0">0</Item>.  <Item name="STATISTICS_MAXIMUM" sample="0">255</Item>.  <Item name="STATISTICS_MEAN" sample="0">108.47642098115</Item>.  <Item name="STATISTICS_STDDEV" sample="0">83.473159026412</Item>.</GDALMetadata>.
GDAL No Data                    : 256
Geo Tiff Version                : 1.1.0
GT Model Type                   : Projected
GT Raster Type                  : Pixel Is Area
GT Citation                     : PCS Name = CF-19 WAC
Geographic Type                 : NAD83
Geog Citation                   : GCS_North_American_1983
Geog Angular Units              : Angular Degree
Geog Semi Major Axis            : 6378137
Geog Inv Flattening             : 298.257222101
Projected CS Type               : User Defined
Projection                      : User Defined
Proj Coord Trans                : Lambert Conf Conic 2SP
Proj Linear Units               : Linear Meter
Proj Std Parallel 1             : 46.6666666666667
Proj Std Parallel 2             : 41.3333333333333
Proj False Origin Long          : -69
Proj False Origin Lat           : 44.41667
Proj False Origin Easting       : 0
Proj False Origin Northing      : 0
Image Size                      : 11888x4135


texts:
I am trying to join points to polygons based on whether their addresses match as well as if the distance between them is less than 2500 feet. I defined a haversine function that computes the distance between the point and the polygon using their coordinates. However I get an error that the sql is invalid.

import arcpy
import harversine
arcpy.env.overwriteOutput=True
arcpy.env.workspace=r"C:\mygdb.gdb"
where="pointshp.CAdd = polyshp.P_Address and harversine.myhaversine(pointshp.NewLong,pointshp.NewLat,polyshp.polyX,polyshp.polyY) < 2500"
keyfield="pointshp.OBJECTID"
arcpy.MakeQueryTable_management(["pointshp","polyshp"],"joined","USE_KEY_FIELDS",keyfield,"",where)
arcpy.CopyFeatures_management("joined","joined")

The error I get:

An invalid SQL statement was used. [joined]

The haversine function I use which is saved in a separate script called harversine.py i.e. it is not a spelling mistake.

import math
def myhaversine(lon1, lat1, lon2, lat2):
    # convert decimal degrees to radians
    lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
    c = 2 * math.asin(math.sqrt(a))
    r = 20887680 # Radius of earth in kilometers. Use 3956 for miles, 20,887,680 feet
    return math.ceil(c * r)


texts: With QGIS adjacent polygons, is it possible to "copy" one polygon edge to a new one? Like a "snap to" or "smooth to" the adjacent polygon?

texts: I'm drawing adjacent polygons in QGIS.  Whenever I create a new adjacent polygon, or if I've moved a polygon that has an adjacent polygon, I need to adjust the borders of the adjacent polygon to be right on the first one.  When adjusting one adjacent polygon I've been using the Vertex Tool, clicking a vertex of the adjacent polygon, and placing it where the first one is or has moved to, one vertex at a time, zooming in for highest accuracy.  Seems perhaps there might be an easier way to do it?  Would a paid product like ESRI have a more elegant solution to this?  I've been a consumer of GIS data, just recently attempting to edit/create.


texts: Joining all hydrographic vectors touching each other (inside same feature) using QGIS?

texts:
I have a complex hydrographic layer where all the hydrographic lines are separate even if it's from the same river.

How do I select all the vectors who are touching themselves without select vector from another river using automation ?

I'm on Qgis.


texts: reprojecting geotiffs to spherical mercator using gdalwarp

texts:
I would like to reproject raster images (tiffs+proj+tfw) files too "google's" projection using gdalwarp.  From what I understand, google uses epsg:900913/epsg:3785.
What I don't know is how to correctly state the source projection.
I believe that UTM32N=epsg:32632.  If I reproject using gdalwarp -s_srs epsg:32632 -t_srs:epsg:900913, would that be correct? If I open that file in global mapper, how do I adjust project and worldfile then?

This is the content of my .proj file:

PROJCS["AUT-UTM32N",
    GEOGCS["GCS_WGS_1984",
        DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["false_easting",500000.0],
        PARAMETER["false_northing",0.0],
        PARAMETER["central_meridian",9.0],
        PARAMETER["scale_factor",0.9996],
        PARAMETER["latitude_of_origin",0.0],

This is the .tfw file:
    2.5000000000
    0.0000000000
    0.0000000000
    -2.5000000000
    550001.2500000000
    5199998.7500000000
    UNIT["Meter",1.0]]


texts: Reproject MOD13Q1 in GDAL using the gringpoints

texts:
I am working with the MODIS MOD13Q1 vegetation product and the reprojection using gdal is fine, however I wish to be more precise and use the GRING to accurately warp the data into wgs 84(epsg:4326). Whilst trying to accomplish this using the resource provided by PolyGeo and Valgur, I tried my best to reproduce the same parameters for my dataset and the results are not satisfactory at all.
For reference, My dataset looks like this

After implementing this function:
def repro_gcp(test_ds,gcp_list):
    #GCP related Informatiom
    gcp_srs=osr.SpatialReference()
    gcp_srs.ImportFromEPSG(4326)
    gcp_wkt=gcp_srs.ExportToWkt()
    test_ds.SetGCPs(gcp_list,gcp_wkt)
    #transformed dataset related information
    dst_srs=osr.SpatialReference()
    dst_srs.ImportFromEPSG(4326)
    dst_wkt=dst_srs.ExportToWkt()
    resampling=gdal.GRA_CubicSpline
    #source dataset related information
    src_srs=osr.SpatialReference()
    wkt=test_ds.GetProjection()
    src_srs.ImportFromWkt(wkt)

    repro_temp=gdal.AutoCreateWarpedVRT(test_ds,src_srs.ExportToWkt(),dst_wkt,resampling,0.125)
    dst_xsize = repro_temp.RasterXSize
    dst_ysize = repro_temp.RasterYSize
    dst_gt = repro_temp.GetGeoTransform()
    repro_temp=None
    dst_ds = gdal.GetDriverByName('GTiff').Create('funtest.tif', dst_xsize, dst_ysize, test_ds.RasterCount)
    dst_ds.SetProjection(dst_wkt)
    dst_ds.SetGeoTransform(dst_gt)
    gdal.ReprojectImage(test_ds,
                        dst_ds,
                        None,  # src_wkt : left to default value --> will use the one from source
                        None,  # dst_wkt : left to default value --> will use the one from destination
                        resampling,
                        0,  # WarpMemoryLimit : left to default value
                        0.125,
                        None,  # Progress callback : could be left to None or unspecified for silent progress
                        None)  # Progress callback user data
    dst_ds=None

    #using the GCPs to transfrom the Hdf file
    xyz=[(74.2179978544028,19.9340873467462),(85.1321831541614,19.8968271260537),(92.39172194296,29.9992439210911),(80.5381230425356,30.0384277080736)]
    row_col=[(11149,0),(1149,6234),(0,6234),(0,0)]
    #Assign GCPs
    gcps=[]
    zp=zip(xyz,row_col)
    for (x,y) , (row,col) in zp:
        gcps.append(gdal.GCP(x,y,col,row))

test_file=gdal.Open("HDF4_EOS:EOS_GRID:MOD13Q1.A2000353.h25v06.006.2018222145654.hdf:MODIS_Grid_16DAY_250m_500m_VI:250m 16 days NDVI")
    repro_gcp(test_file,gcps)

My output looks like this

I am using the cubic spline interpolation here, I would like to know where I am going wrong .


texts: Lambert cylindrical projection support in PostGis

texts:
Does PostGis support Lambert cylindrical projection or is there a way to add support for it? I'm probably searching the wrong terms, but google isn't bringing up anything useful. Any help is appreciated.

EDIT: if not Lambert, are there are other cylindrical equal-area projections which can be used in PostGis.


texts: Using IGeoProcessor2 to update Oracle Spatial

texts:
In the below example the input and out are shape files read from folders.Is it possible to read and write Oracle Spatial Geometry tables? I have ArcGIS,ArcSDE and Oracle Spatial 11g.
How to set parameters to create a new Geometry table if it does not exist?
How to append the data if a geometry table already exists ?

IGeoProcessor2 pGeoProcessor = new GeoProcessorClass();
IVariantArray parameters = new VarArrayClass();
parameters.Add(pFlayer);
parameters.Add(OutLayerPath + "\\" + pFlayer.Name + "MBM.shp");
parameters.Add("CIRCLE");
parameters.Add("NONE");
parameters.Add("");
parameters.Add("MBG_FIELDS");
pGeoProcessor.Execute("MinimumBoundingGeometry_management", parameters, null);


texts: looping the mask in python for arcgis 10

texts:
I am working with raster data on a worldwide level with a resolution of 30 arc seconds.
I need to export the data into a table format but given the high resolution you can not run the "Sample" (Spatial Analyst) command for the entire dataset at once.

Therefore, I have created a mask for each country(with The Set Null command). I was trying to write a code in Python, using a loop where each time the mask will become
a different zone, i.e. env.mask = zone1, then zone2 and so on.
However, I suppose that within a python session one cannot keep on changing the mask using a loop. Could you kindly suggest a way loop the mask (so that I can just have one code for all the countries but my program still runs)?


texts: Artifacts in gdalwarp output (from polar stereographic to azimuthal equidistant)

texts:
I'm trying to warp some GeoTIFF images (from NSIDC), from a polar stereographic projection, to an azimuthal equidistant projection using gdalwarp:

gdalwarp -t_srs 'EPSG:102016' \
  -te -7000000 -7000000 7000000 7000000 \
  input.tif output.tif

The result looks fine overall, except for a single horizontal row of pixels in the middle that looks "shifted" inward toward the center.

Is there any way to avoid this?

Edit: Input image (as a gif) for reference


texts: Creating a Data Cube from point measurements using ODC

texts:
I wonder if Open Data Cube software can also handle point measurements (none-raster information). I am interested to create a data cube for CYGNSS satellite data. Each CYGNSS data file (in netCDF format) include daily point measurements (more than 8500 samples) from different locations. https://podaac.jpl.nasa.gov/dataset/CYGNSS_L1_V2.1
Is it possible to use ODC to resample CYGNSS data to a regular grid (netCDF format optimized in time direction)?


texts: Converting SHP file of river drainage to TXT using ArcGIS for Desktop?

texts:
I need to convert a shp file of river drainage into *.txt to run in another software.

I am in trouble since the SAMPLES in ArcMap toolbox is deprecated after version 10.


texts: Raster spplot legend: how to add a title beside it

texts:
After studying for a while how to create a chart from a raster object, I got almost everything I needed except to add a title to spplot's legend.

What I tried:

to use mtext("XXX (m)", side=4),
to add a sp.text outside the grid,
to use legend.args = list(title="XXX (m)").

The below picture pinpoint, the feature I added manually to my plot, and in which I would like to include using R.

Here is one example, of what I got so far:

require(raster)
require(sp)
require(lattice)

north <- list("SpatialPolygonsRescale", layout.north.arrow(type=1), offset = c(0.95,0.85), scale=0.1)
scale <- list("SpatialPolygonsRescale", layout.scale.bar(),
      offset = c(0.55, 0.03), scale = 0.4, fill = c("transparent","black"))
txt1 <- list("sp.text", c(0.55, 0.08), "0")
txt2 <- list("sp.text", c(0.75, 0.08), "0.2")
txt3 <- list("sp.text", c(0.95, 0.08), "0.4")
raster_layout <- list(north,scale,txt1,txt2,txt3)
cuts <- c(110,120,130,140,150,160,170,180,190,200)
spplot(raster(volcano), scales = list(draw=T), at=cuts, col.regions = palette(gray(seq(0,0.9,len=9))), sp.layout= raster_layout, add=T)

Thank's for your help.


texts: Selecting features in order using python in QGIS

texts:
I am building a plugin in QGIS 2.14 that takes a selection of lines and creates one followong the path of those lines already selected.

My issue is that the traditional selection does not respect the order we selected the features.
I saw  one topic talking on this here : Get features in order they are selected using QGIS but I didn't manage to implement it.
Do you guys have an example of a code using it, as the need is :

Select the features from the map
Click on the plugin button to peform the operation


texts: Using python libary simplejson and urllib in PostGIS

texts:
I installed PostGreSQL9.3 and want to create a Pl/Python function using simplejson library and urllib. First I use OSGeo4W Shell to test the similar python script, it works. Then in PostgreSQL, it doesn't work; the error info is "can't use a string pattern on a bytes-like object".

CREATE OR REPLACE FUNCTION scratch.GetWeather(lon float, lat float)
RETURNS float
AS $$
    from urllib.request import urlopen
    import json as simplejson
    data = urlopen('http://api.openweathermap.org/data/2.1/find/station?lat=%s&lon=%s&cnt=1' % (lat, lon))
    js_data = simplejson.load(data)
    if js_data['cod'] == '200':
        if int(js_data['cnt'])>0:
            station = js_data['list'][0]
            print ('Data from weather station %s' % station['name'])
            if 'main' in station:
                if 'temp' in station['main']:
                    temperature = station['main']['temp'] - 273.15
                else:
                    temperature = None
            else:
                temperature = None
        else:
            temperature = None
        return temperature

    $$ LANGUAGE plpython3u;

Error message is

"TypeError: can't use a string pattern on a bytes-like object
CONTEXT:  Traceback (most recent call last):
  PL/Python function "getweather", line 5, in <module>
    js_data = simplejson.load(data)..."

How to fix this? Thanks.


texts: Link to external object in Google MyMaps

texts:
I have got a map created in Google MyMaps, which contain a lot of points:

https://www.google.com/maps/d/edit?mid=1yyHAfiqoMCB2X3oswH2HKemh304&ll=-2.3807679248749665%2C0&z=2

Is it possible to input a link to every single point in the description, when click "Edit" option?

AN only option, which I can see there is the photo. I would like to have also a link to particular URL from that place.

I was trying to input some HTML stuff such this:

a href="https://www.britannica.com/science/rift-valley"Spectacular Great Rift Valley> but it doesn't work. It looks like is not supported by Google.

An only option is pick up the point at arleady existing location provided by Google Maps, which arleady has some link.

Is the Google MyMaps tool support a link option at all?


texts: In QGIS map composer, how do I get the Item Properties dialog back?

texts:
Yes, well if there is button there, one is bound to click it....
I closed both the Item Properties dialogs in Map Composer.

And I can't get them to give me any joy no more...

QGIS Trunk


texts: QGIS 2.12.1 crashes on Mac OSX-11 because of Incompatible PROJ library version

texts:
I have just installed QGIS 2.12.1 on my mac OSX-11 El Capitan and all relevant libraries - on running QGIS and loading a previous project the program shuts down.
The error reported is:

Dyld Error Message:
  Library not loaded: /Library/Frameworks/PROJ.framework/Versions/4/PROJ
  Referenced from: /Applications/QGIS.app/Contents/MacOS/QGIS
  Reason: Incompatible library version: QGIS requires version 11.0.0 or later, but PROJ provides version 10.0.0


texts: OpenLayers ZoomToExtent doesn´t seem to work properly with WFS layer

texts:
My website has a list of routes and a map. When I click on each route it appears its associated information and the WFS layer is selected. As a base layer I have all the routes obtained with WMS.
This is a portion of the html code:

 <a href="#" onclick="javascript: displayRouteById(245, 417000, 4593700, 426200, 4597475); return false;">

What I can´t get is the zoom to each route although the extents are correct. Here is my OpenLayers code:

var map, route, extent, routeWMS, draw, modify;

function init() {
    var bounds = new OpenLayers.Bounds(417000,4589000,428000,4592000);
    var resolutions = [1100,550,275,100,50,25,10,5,2,1,0.5,0.25];
    var genericMapOptions = {projection: "EPSG:23031", units: 'm', maxExtent: bounds}

    routeWMS = new OpenLayers.Layer.WMS();

    route = new OpenLayers.Layer.Vector(
        "Base Vector Layer",
        {
            isBaseLayer: false,
            strategies: [new OpenLayers.Strategy.BBOX()],
            projection: new OpenLayers.Projection("EPSG:23031"),
            protocol: new OpenLayers.Protocol.WFS({
                version: "1.1.0",
                srsName: "EPSG:23031",
                params: {
                    map: "c:/inetpub/wwwroot/mapServer/maps/stcugat/map.map",
                    USERID: userId,
                    PERFILID: profileId
                },
                url: "http://localhost/mapserver/mapserv.exe",
                featureType: "routes",
                featurePrefix: "ms",
                geometryName: "the_geom"
            })
        }
    );

    var layers = [routeWMS, route];

    route.events.register(
        'loadend',
        route,
        function() {
            // I get the same extents with the three options
            alert(this.getDataExtent() + "\n" + extent + "\n" + this.features[0].geometry.getBounds());
            //map.zoomToExtent(this.getDataExtent());
            //map.zoomToExtent(extent);
            map.zoomToExtent(this.features[0].geometry.getBounds());
        }
    );

    var mapControls = [
            new OpenLayers.Control.Navigation()
    ];

    var mapOptions = OpenLayers.Util.extend({
        resolutions: resolutions,
        controls: mapControls
        }, genericMapOptions);

    map = new OpenLayers.Map('map', mapOptions);
    map.addLayers(layers);

    map.zoomToMaxExtent();
}

function displayRouteById(id, minX, minY, maxX, maxY) {
    extent = new OpenLayers.Bounds(minX, minY, maxX, maxY);
    route.protocol.params.PERFILID = id;
    route.refresh({force: true});
}

Any advice?


texts:
I tried to perform 'Solve location allocation' in ArcGIS and I received an error when I provided input road data to the to Network dataset section. It returned the following error: Error 030033 : Parameter does not contain a network dataset data element

How can I solve this problem?


texts: How to generate schemas of GIS hardware / software / workflows?

texts:
Quite often I find myself in a situation where I need to explain hardware or software GIS infrastructure to people with no or limited knowledge of GIS. Or give an overview of geoprocessing tasks involved in an analysis.

In such cases picture can be trully worth thousand words. Doing last minute preparations I usually relied on quick and dirty PowerPoint slide with random assortment of free icons and images from around the web but I would be more happy to invest little bit more time into something that will look better.

I'm thinking about schemas resembling for instance images like this:

Or this:

What kind of software (preferably free) would be best suited to create such images?

Is there any place where you could get (preferably for free) 'GIS-specific' icons used in such diagrams?


texts:
I have QGIS linked to PostgreSQL but when I add layers from the db onto the project, people can easily see the user by simply mousing over the layer (which is the master user) so I need to hide it.
How to do this?


texts: Creating random points in ArcGIS on polyline in polygon?

texts:
I am trying to create 100 random points, 1 within 500m of 100 points I already have.  I have created a buffer and can easily create 1 random point in each, but now I want those points to lie on a polyline e.g. road in that area.

I have clipped the line layer and can get 1 point per line (road) but I only want one point in the buffer.  If a buffer contains 10 roads then I get 10 random points in that area.  I tried dissolving the roads but then I only get 1 point for all the buffered areas.


texts: Show differenet SVG icons on one layer at runtime

texts:
I need to show several SVG icons (more than 388 icons) with different name on a single layer(if possible on memory layer). I am using QGis 2.8 and qt4.8.6 The icons represent different kind of entities such as: plane, radar, weapons and many more... They will be added to the layer at run time based on some preconditions. I've tried using categorized symboling but it only works for different symbols for different values, correct me if I am wrong. This is what I've tried so far and also tried categorized symboling but can't display the icons the way i want them to be displayed.

  `QString IcoType1 = "Path to my icons";
    QString src = "D:/NYC_MUSEUMS_GEO/NYC_MUSEUMS_GEO.shp";
    QgsVectorLayer *layer2 = new QgsVectorLayer(src,"Symbol Layer","ogr");

    QStringList MyPropertyName;
    MyPropertyName << "fill" << "name" << "outline" << "outline-width" << "size";
    // list of corresponding value 1
    QStringList MyPropertyValue;
    MyPropertyValue << "#f3e54d" << lBoyIcoType1 << "#000000" << "0" << "20";
    QgsStringMap PropertyMap;
    for (int i =0; i < MyBoyPropertyName.size(); i++){
        PropertyMap.insert(MyPropertyName.at(i),MyPropertyValue.at(i));
    }

    QgsSymbolLayerV2* MySymbolLayer = QgsSvgMarkerSymbolLayerV2::create(PropertyMap);

    layer2->rendererV2()->symbols(QgsRenderContext())[0]->changeSymbolLayer(0,MySymbolLayer);
    layer2->rendererV2()->symbols(QgsRenderContext())[0]->appendSymbolLayer(MySymbolLayer);
    QgsMapLayerRegistry::instance()->addMapLayer(layer2,TRUE);
    QgsMapCanvas *canvas = new QgsMapCanvas(this);
    QList<QgsMapCanvasLayer> layers;
    layers.append(QgsMapCanvasLayer(layer2,TRUE));
    canvas->setLayerSet(layers);
    canvas->setExtent(layer2->extent());
    showMaximized();
    setCentralWidget(canvas);
`


texts: Create an animation that rotates around an central axis

texts:
I am looking to create an animation of a 3D model.  I want the model to be the central object in the animation, so rotating around the model, always looking at the model.  I am wondering what is the best way to perform this.  This seems to be a basic operation, but it is becoming more tricky than I anticipated.  I would like to use ArcScene or Google Earth, but I am open to suggestions.

Thanks...


texts: Python, comtypes and ArcObjects: Error creating AppROT object

texts:
I am experimenting with comtypes and ArcObjects under Python 2.6.5 and ArcGIS 10 SP1. I'm using the pure Python method to wrap the ArcObjects OLBs described in this answer, but getting an error in the comtypes.CoCreateInstance method.

Here is the code I am running:

def WrapModules():
    #force wrapping of all ArcObjects libraries (OLBs)
    import os
    import comtypes.client
    # change com_dir to whatever it is for you
    com_dir = r'C:\Program Files\ArcGIS\Desktop10.0\com'
    coms = [os.path.join(com_dir, x) for x in os.listdir(com_dir) if os.path.splitext(x)[1].upper() == '.OLB']
    map(comtypes.client.GetModule, coms)

def GetApp():
    """Get a hook into the current session of ArcMap"""
    from comtypes.gen import esriFramework
    pAppROT = NewObj(esriFramework.AppROT, esriFramework.IAppROT)
    if pAppROT is not None:
        iCount = pAppROT.Count
        if iCount == 0:
            print 'No ArcGIS application currently running.  Terminating ...'
            return None
        for i in range(iCount):
            pApp = pAppROT.Item(i)  #returns IApplication on AppRef
            if pApp.Name == 'ArcMap':
                print "ArcMap found"
                return pApp
        print 'No ArcMap session is running at this time.'
    print "No AppROT found"
    return None

def NewObj(MyClass, MyInterface):
    """Creates a new comtypes POINTER object where\n\
    MyClass is the class to be instantiated,\n\
    MyInterface is the interface to be assigned"""
    from comtypes.client import CreateObject
    import traceback
    try:
        ptr = CreateObject(MyClass, interface=MyInterface)
        return ptr
    except:
        print traceback.format_exc()
        return None

if __name__ == "__main__":
    WrapModules()
    pApp = GetApp()
    if pApp is not None:
        print "HWND: %d" % pApp.hWnd
    else:
        print "No ArcGIS application found!"

And here is the output from the script:

Traceback (most recent call last):
  File "C:\temp\ComHelpers.py", line 35, in NewObj
    ptr = CreateObject(MyClass, interface=MyInterface)
  File "C:\Python26\ArcGIS10.0\lib\site-packages\comtypes\client\__init__.py", line 235, in CreateObject
    obj = comtypes.CoCreateInstance(clsid, clsctx=clsctx, interface=interface)
  File "C:\Python26\ArcGIS10.0\lib\site-packages\comtypes\__init__.py", line 1145, in CoCreateInstance
    _ole32.CoCreateInstance(byref(clsid), punkouter, clsctx, byref(iid), byref(p))
  File "_ctypes/callproc.c", line 925, in GetResult
WindowsError: [Error -2147221231] ClassFactory cannot supply requested class

No AppROT found
No ArcGIS application found!

Thanks for any insights you might have!


texts: Leaflet MarkerCluster plugin. Weighted Points

texts:
I have an huge amount of points I'd like to have on the map so instead of sending millions of points across I have grouped together points that are "close enough" to be the same point. So I have moved from 1 000 000 + points down to 10 000s of points which is working well.
Now however I now no longer have an accurate representation.

So ideally I'd like to a way to modify the child count. I was thinking about overriding the create icon class and looping through all the individual icons getting the weighting for each point and tallying it up. Then at the end have

cluster.getChildCount() + total_weighting

but that seemed inefficient and was wondering if anyone new a better approach.


texts: Why there is a shif in MODIS coordinate (compared with Landsat/google earth image) after mosaicing and subsetiing the study area?

texts:
I am facing some difficulties in processing MODIS data set.

To have MODIS data over my study area, I need to mosaic two tiles of MODIS (h29 v12 and h30 v12). I have used the MRT tool by NASA (the website based) to reproject the MODIS data from sinusoidal to UTM (zone 55s). And used VISAT (software by ESA) to mosaic the data.

It is not clear for me why the final subset over the study area comes with some coordinate shifts (13') when it is compared against google earth and land use map derived from Landsat imagery. I guess the mosaicing process is the main reason for this coordinate shift, but I heard that it could be because of reprojection process.

Can you give me some clue about what could be the reason for this coordinate discrepancy or how to solve this problem?


texts: Iterate over elements in Rasterio window and obtain coordinates using transform.xy

texts:
Does a window generated by rasterio maintain the dataset origin and absolute col/row indices?
I'd like to be able to iterate over a window of a dem dataset, find the highest value in the window and return it's lat long coordinates using rasterio.transform.xy or similar.
Obviously the below isn't suitable (the arr is invalid), it illustrates what I'm trying to achieve. Window is generated using slices and row/col high and lows which are provided by a valid function called extract_block.
with rasterio.open(inputraster) as dataset:
        for (x,y) in zip(*nonzero):
            row_low, row_high, col_low, col_high = extract_block(arr,[x,y],windowsize)
            rasteriowindowarr = dataset.read(1,window=Window.from_slices(slice(row_low, row_high), slice(col_low, col_high)))
            with rasterio.open(rasteriowindowarr) as windowarr:
                for i,j in np.ndindex(rasteriowindowarr.shape):
                    coordinatex,coordinatey = rasterio.transform.xy(windowarr.transform, i,j,offset='center')
                    print ("long of result is " + str(coordinatex))
                    print ("lat of result is " + str(coordinatey))


texts: Cannot get Remove From Current Selection to work in ArcMap?

texts:
I am a Newbie using ArcGIS 10.3 for Desktop:

After importing a whole watershed streams Layer, I selected by location a subset of the streams which was already defined with a separate name. "intersect the layer feature", thus the subset streams become highlighted
To take out a piece of the selected subset I did the following:

got closet to the area to be extracted using the + zoom tool
activated the "select features" button for an interactive selection
ensured that only one layer become affected by the changes:
Table of Contents, Make This The Only Selectable Layer,
-Then: Selection, Make This the Only Selectable Layer
Then: Selection, Remove From Current Selection
clicked and dragged small box over a part of the feature... and there was no extraction.
What did happen was that the subset lost its highlighted status which depicted the streams that were related with the larger watershed.


texts: Make a CSV/spreadsheet showing number of polygons in polygons?

texts:
I have a QGIS project containing two multipolygon layers: 2010 Census block groups and Parcels. Each layer is in a separate shapefile.

The thick-lined layer represents block groups. The thin-lined layer is parcels. The block group layer attribute table shows that each block group has a unique ID.

How do I output the number of parcels in each block group into a spreadsheet-style format? Is it possible in QGIS or must I choose another program?


texts: Joining fields which must still be created in ArcGIS ModelBuilder?

texts:
I am in the process of creating a generic model in ArcGIS Model Builder. At one point in my model, I have to do a join between an existing layer (entry parameter) and a table issued from the "Zonal Statistics as table" tool.

I get an error when I launch the model from the Toolbox, as the field to do the join doesn't exist. When I create a "Zone field", "Input join field" and "Output join fields" as parameters, I can get the fields to be defined by the user. However, he has to write the name of the "Zone field".

Is it possible to select that value in a list (like the fields from my existing layer)?

If I am doing things wrong, how can I get the model to say which fields to do the join on?

One layer is a .shp that is defined as a parameter. The other is a layer produced by a raster calculator in the geodatabase then reclassed and to which a "Zonal Statistics as table" was applied.

Version 1 that returns an error:

Version 2 that works, but doesn`t enable the user to pick the field the join needs to be on in a list:


texts: Serving large (1.7Tb) dataset of aerial images with GeoServer

texts:
I have a large dataset of aerial images consisting of 8000 TIFFs where single images are no larger than 0.5G and all have have .ovr pyramid files with them.

Previously I have served this dataset as imageMosaic, but I'm wondering should I use imagePyramid instead. At what point is imagePyramid better than mosaic of TIFF+ovr or is it ever?

As far as I understand GeoServer imageMosaic supports ovr-files so is there any benefit of creating pyramid folder structure?


texts: Determining if two polygons are overlapping in Google Maps?

texts:
I'm working with google maps and polygons, but I have to validate that polygons are not overlapping each other.

I have a function that returns TRUE when a point (lat, lng) is inside a polygon, but it's not enough to determine if every point of a polygon is inside another polygon.

Any suggestion?


texts: Joining attributes in 100s SHPs back to main SHP file batch in QGIS

texts: I have ~200 shapefiles with data in neighborhoods by census blocks across the city. I want to join each of these shapefiles back to the main census block file so that a count column in each of the neighborhood shapefiles gets joined back to the original city census blocks shapefile. I thought I could do this in a batch process using the "join attributes table" function but it seems like this batch process will create 200 joined files, instead of just appending each joins onto the city shapefile table. Does anyone know how I can do this in QGIS? I am using QGIS 2.14.3


texts: How to apply symbol in group to map attributes in qgis 2 using osm roads as an example

texts:
Update: While this question does point to available style guides my question is more how to go about similar operation by yourself and associating your own symbols with osm maps

I downloaded OSM based maps from GeoFabrik and my question is how to associate the symbols shown in the screen shot with the fclass attribute automatically so vectors of ftype type "Motorway","Primary" etc all use the appropriate symbol
I tried changing the symbol type to Categorized by ftype but that still meant I had to manually edit each one and pick the correct symbol.

p.s. are the symbols included in the osm file or are they common to qgis?


texts:
i need world hillshade with resolution about 1x1 km.

Do you know where i can it download ?


texts: Merging features which share geometry but not features

texts:
I'm trying to create one vector layer out of two, one of which overlies the other.  I've merged the two layers but now I just have a load more polygons which share no attributes.  How do I merge the features which have the same (or overlapping) geometry while retaining the attributes of both?
Background - I want to add 'ancient woodland' and all its attributes to 'all woodland' and all its attributes.  Clearly 'ancient woodland' is a part of 'all woodland' so it always shares geometries.


texts: ArcMap, ArcCatalog slow to open on laptop with ample resources?

texts:
I have a recent installation or ArcGIS Desktop on a 64 bit Windows 7 Machine, with 8GB ram, intel core i7 processor, and a solid state drive.

It takes ~2 min to open either ArcMap or ArcCatalog from the start menu (e.g. without loading any data).

As a new installation, I have only the default folder connections.

Is this normal, or is it possible that there is something wrong with my installation?


texts: OGR error creating geometry for feature when deleting nodes

texts:
I have a line Shapefile that I need to remove many nodes from.

Some I can delete fine, others delete but when I hit save I receive the following error:
OGR error creating geometry for feature

Using QGIS version 2.14.0-Essen in Win10 64 bit.

Any ideas? I've converted it to multipoint then back to single point with no change.


texts: Bing Map not rendering in New Print Layout QGIS 3.0

texts:
I am trying to sort out a map with a Bing aerial map plugin, yet when I select New Print Layout, the aerial image will not load/render, yet all other layers do.

Can anyone point me in the right direction of how I can sort this issue out?

I am pretty new to QGIS...


texts:
Using javascript in GeospatialPDF, is it possible to zoom in to an extent, say when the document is opened the first time ?

To what extent can I use javascript to access the geospatial data? pan, zoom, turn-on/off layer ?

From GDAL Geospatial format page:
Create of PDF with some JavaScript :

gdal_translate -of PDF my.tif my.pdf -co JAVASCRIPT_FILE=script.js

where script.js is :

button = app.alert({cMsg: 'This file was generated by GDAL. Do you want to visit its
website ?', cTitle: 'Question', nIcon:2, nType:2});
if (button == 4) app.launchURL('http://gdal.org/');


texts: Split layers by field QGIS version 3.2.3

texts:
I want to split a vector layer based on values in a specific field.
In older versions I could use the plugins 'LayersbyField' or 'Split Vector Layer', but both seem not to be available anymore in version 3.2 of QGIS.

Are there quick and easy alternatives?
I don't like to 'select' and 'save as' for every value, since that takes too much time.


texts: Are objects on periodically updated satellite imagery expected to shift and what to do about it?

texts:
Background 1:
The organization I work at creates and maintains the road network of the entire country as line features, along with some points of interest. They do this by acquiring high resolution color imagery for the country and digitizing over it. When a new updated satellite imagery is acquired (from Airbus or Digital Globe), the digitized features look shifted.

Background 2:
Being a frequent use of the Waze app, I sometimes use the waze editing tool to update streets, landmarks, parking lots etc. The Waze editing tool provides google maps to digitize over and update features on the app. When google periodically updates its satellite imagery, all the digitized streets and features are shifted. This shift doesn't affect the usability of the waze app, but when using the waze editing tool, can seem to indicate the digitization was done incorrectly.

The two part question I have is:

1) Are two satellite images of the same area always expected to have a shift because of the way satellite imagery is captured, processed and provided by imagery providers?

2) Is there any solution that the imagery customer can take or the imagery provider can take to ensure every imagery of a given area does not show a shift?


texts: arcsde huge size file upload into Oracle spatial

texts:
I have 35million records (tomtom multinet) of file gdb file in my desktop machine. I need to upload this file into Oracle spatial as SDO_geometry file.

I already tested the connection with small size file gdb file  and it is success.

Is there any best practice here to upload this huge file (from a desktop machine -windows 7 to oracle ?

my env: arcgis desktop 9.3.1 (sp2) and arcsde (arccatalog and arsde are in same machine), oracle 10g (db is in different location).

Or, whether i have to upload it as splittted files (5m rows each) into oracle spatial, then doing append in a file through arccatalog. I have default dbtune.sde file, commit interval 1000.

Pls. share your experience in this.


texts:
I created a Python plugin in QGIS and I want to make it public, so that any user could download and use it

I serached in https://plugins.qgis.org/plugins/ , but to upload the plugin , I have to have an OSGEO acount. I looked in https://id.osgeo.org/ldap/create, apparentelly I can not create it the acount as usual, I sent the request to mantra-request@osgeo.org but it is taking too much time.

Is there any other way to upload my plugin ?


texts:
GeoServer version 2.18.0 running under Java 1.8.0-131 on a Windows machine.
I am trying to use the GeoServer REST API to seed GeoWebCache. I can get the seed process to commence and can follow its progress using the commands documented. My problem is that I am unable to stop the seeding process once it has started. I have tried the following methods as detailed in the seeding and truncating section of the manual.

Kill the task for a specific layer. POST command of the form :

http://localhost:8080/geoserver/gwc/rest/seed/MYStore:MYLayer
The above line returns 400 - bad request
Various attempts to provide some body for the POST also fail. I have tried JSON of the form
{ seedRequest : { "name" : "MYStore:MYLayer" }}  or the equivalent XML.

Kill all tasks

http://localhost:8080/geoserver/gwc/rest/seed
This returns - Chunk[] is not a valid entry

Using the kill all option as documented (curl attempt shown)

curl -v -u username:password -d "kill_all=all" "http://localhost:8080/geoserver/gwc/rest/seed"

This returns chunk [] is not a valid entry
I have no difficulty kill the seeding task using the GeoServer web interface.
What am I doing wrong?


texts: Identifying coordinate format from exif data

texts:
I am using the Python exif module to extract coordinates from photos. The methods for doing so are: gps_latitude and gps_longitude.
The output for latitude is: (32.0, 37.0, 37.15) and longitude is: (116.0, 30.0, 37.0)

I extracted the location using another tool as well:

32°37'37.2"N 116°30'37.0"W and 32.6270000, -116.5102778

How can I identify the output format: (32.0, 37.0, 37.15), (116.0, 30.0, 37.0) so I can covert it to perhaps decimal degrees?
from exif import Image

with open('IMG_0252.jpg', 'rb') as image_file:
    my_image = Image(image_file)
    print(my_image.gps_latitude, my_image.gps_longitude)


texts: Cesium: How to retrieve the UUID of the building model onclick?

texts:
I have 3D building models, being loaded via JSON with gltf format. I would like to retrieve the individual model's uuid when the user clicks on the building.How can I do it on cesium platform?
Example of UUID : _amk_blk400b_BD.X7jqjezPgOHx6mfXYwv7.


texts:
Are there hosted nominatim services that can do all it's duties (particularly reverse-geocoding) all over the world?

We have a high volume of lookups (we dont want to be throttled) but we want to avoid the hassle of hosting Nominatim ourselves.


texts: Use update cursor in arcpy to count number of points within each polygon

texts:
I have two shapefiles. One contains many points across a state. The other is all the counties within the state. I would like to do two things:

use a cursor to count how many points are in each county and then populate a new field with those data; and

create a true/false field in which true = the county contains more than two points.

I started the first one but was unsuccessful, so the second I'm not sure how to accomplish. This is what I have so far...

import arcpy

points = r"C:\Users\points"
counties = r"C:\Users\counties"
#create new fields
arcpy.AddField_management(counties, "numOfPoints", "LONG", field_length = 10)
arcpy.AddField_management(counties, "hasTwoPoints", "LONG", field_length = 10)

points_lyr = arcpy.MakeFeatureLayer_management(points,r'in_memory\points_lyr')
counties_lyr = arcpy.MakeFeatureLayer_management(counties,r'in_memory\counties_lyr')
# Create an update cursor to access and update states features
#NAME has names of counties and numPoints will be populated with number of points
fields = ['NAME','numOfPoints','hasTwoPoints']

with arcpy.da.UpdateCursor(counties_lyr,fields) as cur:
    for row in cur:
        county = row[0]
        where = '"NAME" = \'{}\''.format(county)
        arcpy.SelectLayerByAttribute_management(counties_lyr,'NEW_SELECTION',where)
        arcpy.SelectLayerByLocation_management(points_lyr,'INTERSECT',counties_lyr)
        number_of_points = int(arcpy.GetCount_management(counties_lyr).getOutput(0))
        row[1] = number_of_points
        cur.updateRow(row)
print('Operation complete.')

The current script works, but it populates the field 'numOfPoints' with all 1 values which is not correct.

then at the end to add the true/false based on the output to be a stand alone script,

fields2 = ['numOfPoints','hasTwoPoints']

with arcpy.da.UpdateCursor(counties_lyr,fields2) as cur:
   for row in cur:
       if row[0] > 1: row[1] = 'TRUE'
   else: row[1] = 'FALSE'
   cur.updateRow(row)

I get this error: StopIteration: iteration not started


texts: I am a newbie in PostGis and pgRouting. I have installed and configured both on my windows PC, following FOSS4G tutorial. I am able to see the result in browser, using GeoServer and OpenLayers. However, the starting and ending points are not the nearest edge of road, but the nearest node. I have created a 'Source', Target' and 'Edge' filed in my table. Using 'assign_vertex_id' function, source and target fields are getting populated. What to do with edge filed? Am I missing any step?


texts:
I work with QGIS 3.12 and the default pan map symbol is

I try to find a way to change it into an arrow.


texts: version of PostgreSQL/PostGIS supported by Embedded Cadcorp SIS 8.0

texts: Can anyone tell me where I can find out what versions of PostgreSQL/PostGIS are supported by embedded Cadcorp 8.0 which I believe is based on Cadcorp SIS Map Modeller 8.0.


texts: How to handle list variable in a where clause

texts:
I am calling a list of names in from a text file and would like to use each of those to make new layers. I am running into an issue where I get "TypeError: cannot concatenate 'str' and 'list' objects" and am not sure how to get past this. I tried a couple of different where clauses without much luck.

*EDIT: Working code below. The issue was I am calling the file directly in the walk and the list fields. so when I am appending the path I needed to remove the filename from os.path.join(dirpath, filename).

import arcpy, os
from arcpy import env

#setup input varible
ws = r"M:/"

cl = r"/Admin_1.shp"

dem = r"M:/SRTM/dem"

#set output as a variable
out = r"M:/Elevation_Masks/Masks"
    with open('M:/Elevation_Masks/test_space/elevCountries_Test.txt', 'r') as f:
    cb = [line.strip() for line in f]

CNTRY_NAME_FIELD = 'CNTRY_NAME'

for dirpath, dirnames, filenames in arcpy.da.Walk(cl, datatype="FeatureClass", type="ALL"):
    for filename in filenames:
        if CNTRY_NAME in [f.name for f in arcpy.ListFields(cl)]:
            for country_name in cb:
                print "Country:", country_name
                #setup where clause
                where = '"%s" = \'%s\'' % (CNTRY_NAME, country_name)
                arcpy.Select_analysis(os.path.join(dirpath), r"M:/Elevation_Masks/test_space/Country_Masks_Test/clipLayers_Test/clipLayers_Test.gdb/" + country_name + "_Boundary", where)


texts: Converting arcs to straight line segments?

texts:
I have a large number of geometries, currently stored in Oracle, some of which contain arc segments.

I need to convert these geometries into LRS (Linear Referencing System) format so I would like to convert all the arcs to straight line segments.

I may have to write something in PL/SQL but, being a lazy sort, I was rather hoping someone had done this before.

Doesn't have to be in PL/SQL, I also have ArcGIS for Desktop, MapInfo, GeoMedia and FME if any of those would do it.


texts:
I'm trying to pull a LineString geometry out of PostGIS and parse it with OGR (python bindinds).

from osgeo import ogr
import psycopg2

connection = psycopg2.connect("...")
cursor = connection.cursor()

query = "SELECT geom FROM points LIMIT 1"

cursor.execute(query)
row = cursor.fetchone()

wkb = row[0]
geom = ogr.CreateGeometryFromWkb(wkb)

cursor.close()
connection.close()

I've tried:

wkb = bin(int(row[0], 16))

and:

SELECT ST_AsEWKB(geom) FROM points LIMIT 1

OGR does not want to parse it. Keeps giving the following error:

ERROR 3: OGR Error: Unsupported geometry type


texts: GPL License - App only uses QGis components as a GIS display

texts:
I'm still not clear on QGIS licensing ...

I am contracted to write a system that has nothing to do with GIS and GIS analysis, but goes around managing animal tags, and I want to use QGIS to draw on a panel in my app, and display a map of where the animals are as a secondary function.

The app is for a specific client and will not be publicly distributed (sold and installed on client specific computers only and not sold as a program to the public). Does my program still have to be GPL, which will be a problem since there is a lot of other private copyrighted code involved, or do I only have to provide QGIS source code on request?


texts: How to alter this PostGIS ST_distance_sphere query to give the answer for all points in the table, not just one?

texts:
I have the following query which returns the expected result for a single point.  I would like to alter this query to return the answer for every point in the table.

SELECT count(case when(
                   ST_Distance_Sphere(
                       (select the_geom
                        from brc_property where property_r=20), the_geom) < 500000)
                   then 1 else null end)
FROM brc_property


texts: Locating British National Grid Eastings/Northings in ArcGIS Pro i.e. ArcMap Go To XY equivalent?

texts: In ArcGIS Pro is there a way to search for Eastings and Northings for Britain, similar to the Go to XY in ArcMap?


texts:
Based on the answer of sgillies I want to use Shapely's interpolate() to find the points at which to split. But I experience that split() is sensitive to the points being precisely on the geometry you want to split. How to cope with these Floating point issues.

I introduce the following function which should return line segment of equal length:

def get_linesegments(line, n):
    points = MultiPoint([line.interpolate(i/n, normalized=True) for i in range(1, n)])
    print(points.wkt)
    return split(line, points)

Now generate a line:

line = LineString([(0,1), (1,0), (2, 1), (1, 2)])

When the line in split in three parts, the function works correct:

>> lines = get_linesegments(line, 3)
MULTIPOINT (1 0, 2 1)
>> lines.wkt
'GEOMETRYCOLLECTION (LINESTRING (0 1, 1 0), LINESTRING (1 0, 2 1), LINESTRING (2 1, 1 2))'

But when the line is split into 2 parts, we get floating points and the function returns the complete line:

>> lines = get_linesegments(line, 2)
MULTIPOINT (1.5 0.5000000000000001)
>> lines.wkt
'GEOMETRYCOLLECTION (LINESTRING (0 1, 1 0, 2 1, 1 2))'


texts: Zonal statistics as a table error code 010429

texts:
I am trying to extract majority land-use type from a land-use raster (TIFF, unsigned integer, 30X30 m cell size) into a multipart polygon feature class (which contains ~1300 200-m radium circular polygons) using "zonal statistics tool" of ArcGIS 10.3.1.  My goal is to know the major land-use type in each of my circular polygons.

The land-use raster's extent is the conterminous US. I set up the geoprocessing environment correctly (processing extent, ell size), yet I keep getting the following sets of error messages:

.

The detailed error messages are also posted here. I think the issue starts with the first error- there's nothing wrong with the raster dataset and I am not memory full either....

Can anyone help me out?


texts: How to generate a list of unique years from date field?

texts:
I have data that was collected over many years and I'd like to be able to run a python script tool in ArcToolbox (ArcGIS 10.1) to get a list of unique years from the date field.

I found an example of how to format the date field, but this example lists the year from every record, and I just want a list of unique years.

# Imports
import arcpy
from datetime import datetime

# Input Table
fc = arcpy.GetParameterAsText(0)

# Date Field
field = arcpy.GetParameterAsText(1)

rows = arcpy.SearchCursor(fc)
for row in rows:
    datetimeVal = row.getValue(field)
    formattedTime = datetime.strftime(datetimeVal, "%Y")
    arcpy.AddMessage(formattedTime)

I can get a list of unique dates using a search cursor, but I just want the year and not the entire date field:

# imports
import arcpy
from datetime import datetime

# Input Table
inputTable = arcpy.GetParameterAsText(0)

# Date Field
field = arcpy.GetParameterAsText(1)

# Use SearchCursor with list comprehension to return a
#  unique set of values in the specified field
values = [row[0] for row in arcpy.da.SearchCursor(inputTable, (field))]
uniqueValues = set(values)
arcpy.AddMessage(uniqueValues)

How can I get the date in the example above, format it to just get the year, and then run the search cursor to get a list of unique years?


texts: How can I get the field index of the new layer in my script using Python?

texts:
I am trying to use layer.changeAttributeValue(fid, fieldIndex, value) to add values to a field in QGIS

Basically, I'm creating a plugin which will, if a checkbox is ticked, create a new field, then iterate through the dataset, filling that field with the x (easting) coordinate of each point (I also have another checkbox for the same for the y coordinate). As such, I cannot set the field index beforehand because I don't know which will be ticked.

So - how can I get the field index of the new layer in my script using Python?

Also (sorry, I'm new to this), how can I get the fid of the current feature as I iterate through the dataset?


texts: Using ogr2ogr to Export from Access to CSV (or Oracle)

texts:
I received a an Access database too many observations to simply export as a CSV by either the export buttons in Access or copying and pasting.  I tried using ogr2ogr and only get the headers in my "test.csv" file.  What am I doing wrong here?

Conversely, if possible, I'd like to import the table into Oracle.  Is there a method to do that straight from Access?

   ogr2ogr -update -overwrite -nlt none -progress -f CSV 'G:\Office\Heritage\KDFWR SWG Grant Work\2016_data_share\test.csv' 'G:\Office\Heritage\KDFWR SWG Grant Work\2016_data_share\Data Share\ALL_STWG.mdb'


texts: ArcGIS 10.1 Network Analysis- How do I get appropriate network dataset?

texts:
I have acquired London road dataset from the Open Street Map however after creating and building the network dataset in ArcCatalogue the solve tool could not return the routes because of errors no valid stops. That after I had enter the addresses as stops, No result was returned.

Olawumi


texts: QGIS Atlas Grouping Features by Attribute

texts:
I've got the Atlas function down (AWESOME!) and am generating maps for every feature.

What I need to do now is generate maps for groups of features based on a shared attribute.

fid | name | group
01  | a    | 1
02  | b    | 1
03  | c    | 2
04  | d    | 2

Essentially, two maps, one with "zoom to extent" for all group 1 features, and the other zoomed to extent for all group 2 features.

Is possible?


texts: PostGIS TIGER reverse_geocode() returns state only

texts:
I have a clean install of PostGIS 2.1 in Postgres 9.3 on Ubuntu, using apt-get and postgres extensions.

Ran the data loading scripts from Waiting for PostGIS 2.1 - Install PostGIS Tiger Geocoder as an Extension.

normalize_address works well, but geocode returns no result.

Interestingly, reverse_geocode is able to find the state, but no finer detailed information.

  SELECT pprint_addy(r.addy1) As st1, pprint_addy(r.addy[2]) As st2,
  pprint_addy(r.addy[3]) As st3, array_to_string(r.street, ',') As
  cross_str FROM reverse_geocode(ST_GeomFromText('POINT(-71.06941
  42.34225)',4269)) As r;

returns:

  st1   st2 st3 cross_str

  MA

Is this a data loading issue or a software install issue?


texts:
How can I adjust the online basemap imagery file to our existing ArcMap data. Below are steps I perform to add the basemap data.

I Added basemap from "Add Data From GIS Online" in the file menu and
then selected Imagery  and added to the Desktop.
I already have some layer files in our current project in arcmap the details of our

Qatar GCS (Datum)
:

Qatar Projection: Transverse_Mercator

False_Easting: 200000.000000

False_Northing: 300000.000000

Central_Meridian: 51.216667

Scale_Factor: 0.999990

Latitude_Of_Origin: 24.450000

Linear Unit: Meter

GCS_International_1967

Datum: D_International_1967

====================================================================

and the imagery from online data GCS(Datum) :

   WGS_1984_Web_Mercator_Auxiliary_Sphere

   Projection: Mercator_Auxiliary_Sphere

   False_Easting: 0.000000

   False_Northing: 0.000000

   Central_Meridian: 0.000000

   Standard_Parallel_1: 0.000000

   Auxiliary_Sphere_Type: 0.000000

   Linear Unit: Meter

   GCS_WGS_1984

   Datum: D_WGS_1984

Now the problem is the added imagery data is looks like it's shifted around 100m difference with our arcmap hole data.

then I tried transformations to adjust with our arcmap data by selecting the same as our Datum: D_International_1967 in the dataframe properties but not changed anymore. and then i tried different transformations named as Qatar but same result not aligned the data.

am I trying to do the transformations at the right place.?? and any thing i missed here.???

**Right click in the arcmap\Datafareme properties\Coordinate system tab\Layers\selecting the desired layer to transform\transormations\convert from (selecting existing GCS)\into (Selecting desired GSC)
hit Ok.**


texts: How to show on Leaflet: WMS on top of OpenMapTiles (vector PBF)?

texts:
I made a simple Leaflet page with the following content:

Overlays:

a geojson polygon file
WMS polygon layer from Geoserver of my own (tiled=true)

Baselayers:

All 4 Google Maps: road, terrain, satellite, hybrid (from internet); with Google Mutant Leaflet plugin
OpenStreetMaps (from internet)
OpenMapTiles (of my own from my intranet)
a blank baselayer

Everything is smooth but one thing: every time i turn on my OpenMapTiles (OMT), my WMS dissapear from map; and the WMS will be shown normally again with baselayer other than OMT.

I tried to remove Google baselayers (as well as script reference to Google Maps API accordingly) from the page, hence only OMT, OSM and blank baselayer remain. But the WMS won't show with OMT baselayer.

Any clue why this is happening and how to fix this?

The stack includes:

Leaflet 1.3.3
GeoServer 2.7.0 (for the polygon WMS)
Tomcat 8.5.32 on port 8080
OpenMapTiles on IIS behind an F5 reverse proxy


texts:
I am using QGIS 2.18.15. I got some trouble to add an asc-file as a raster layer into QGIS. I tried to add the raster layer the normal way via the tool bar bottom. It did not work. (“The file source is not supported”)

Then I tried to convert the data via Raster: Conversion: Translate (convert format) menu option as it is suggested at StackExchange (Manipulate .asc files in qgis). It did no work ether. Here is the same problem “not recognized as a supported file format.”

The data I would like to add is from the following link. ftp://ftp-cdc.dwd.de/pub/CDC/grids_germany/multi_annual/solar/

It does not matter which data package I downloaded neither of them worked out.


texts: Is it possible to lock records or attributes indefintely if a field is populated?

texts:
I am creating a c# ArcGIS addin where I have a feature class stored in Oracle, when a feature is created some of the fields within the table will be populated. If field X is populated, then certain fields will be locked, if field Y is populated then some other fields will be locked and if field Z is populated then the geometry will be locked.

Locked means no  edits can be done (geometry wise or the attributes) either  within the addin and outside in ArcGIS or even in Oracle.

Is this something that can be done? If so, does anyone have any pointers to get me started?

EDIT - There may be a scenario where the feature class is joined to table and that table will have a field that locks the geometry. How would it be possible to lock the spatial table?


texts: Getting elevation point from DTED file using ArcGIS Runtime SDK for WPF?

texts:
I am still new to the ArcGIS Runtime SDK for WPF.

My project requires me to get the elevation point offline (local) from a DTED file. I had raster map view using mapview, then I want to get the elevation point when I point the mouse tip on raster map that I view.

This is my code so far:

string parentFolder = @"C:\Data\DTED0";
List<string> files = Directory.EnumerateFiles(parentFolder, "*.dt0", SearchOption.AllDirectories).ToList();

Scene.BaseSurface = new Surface();
RasterElevationSource rasterElevationSource = new RasterElevationSource(files);

Scene.BaseSurface.ElevationSources.Add(rasterElevationSource);

Do I need to create new user control for sceneview or just use the that i read raster map on mapview?


texts: Google Places Autocomplete Limit Suggestions to Local Area

texts:
I'm not sure if this is the right place to ask, but can I limit the suggestions that Google Places Autocomplete find to say 5 km of a certain location? If I'm in Dublin I don't want to see suggestions in Canada!

https://google-developers.appspot.com/maps/documentation/javascript/examples/places-autocomplete


texts: Why are the Google basemaps no longer appearing in QGIS?

texts:
I had done a few projects using Google Physical/Streets as basemaps but now when I re-open those projects all layers appear except for the google basemap (other layers include .kml).

I had previously had this issue when I had a poor internet connection - the basemap would not load - but I now have a good connection so that should not be the problem.

I have also tried starting a new project and adding Google Physical there but likewise nothing appears. In the old projects I have tried removing and re-adding the basemap but no luck.

I have tried uninstalling and re-installing the OpenLayers Plugin (version 1.3.6) but the problem persists. When I use the plugin to add OSM however it does work (but I need Google for these projects).

I have verified the CRS and have it set to WGS 84/Peudo Mercator EPSG:3857

I am using QGIS 2.8.2-Wien on Mac

How do I make Google Physical visible again?


texts: VBA Labeling works on some symbol classes but not other

texts:
VBA Labeling formatting works on some symbol classes but not others.

For some symbol classes it will work but others it won't register the formatting and I can't figure out why!

CLR red = '56' green = '168' blue ='0'>" & [LOCA_ID] & /CLR

This is the section of VBA not working, obviously it is closed off in the VBA statement.

[LOCA_ID] is string.

Updated:

Function FindLabel ([LOCA_ID],[SPEC_DPTH],[ERES_NAME],[ERES_RTXT],[ERES_RUNI])
  FindLabel = "<CLR red ='255'>" & [LOCA_ID] & "</CLR>" + " / " + "Sample Depth - " + [SPEC_DPTH] & vbnewline & [ERES_NAME] + " - " + [ERES_RTXT] + [ERES_RUNI]
End Function

[

Any ideas??


texts: Get a LonLat point by distance or between 2 lonlat points

texts:
I know two latlng points and also can get distance between 2 points. But how do I can get a latlng point on 1/3 of the distance of 2 points.
latlng1  = [37.7052918, -8.8420179]
latlng2  = [37.7062884, -8.8425087]

var point1 = new OpenLayers.Geometry.Point(latlng1.lon, latlng1.lat);
var point2 = new OpenLayers.Geometry.Point(latlng2.lon, latlng2.lat);

var line = new OpenLayers.Geometry.LineString([p1, p2]);
alert(line.getGeodesicLength(new OpenLayers.Projection("EPSG:4326")));


texts: How do I keep my python loop going if there is an error with an SDE workspace?

texts:
I have a folder with several .SDE connections.  I want my Python code to execute various SDE maintenance procedures (such as compress) on each of these databases.  If there is a bad/corrupt/not permissioned .SDE file in the workspace, the process doesn't execute on any of the other databases either. How do I get it to attempt to compress each database in the workspace, even if one fails?

import arcpy

arcpy.env.workspace = "F:/Scripts/DatabaseConnections"

workspaces = arcpy.ListWorkspaces("*", "SDE")

try:
    for workspace in workspaces:
        arcpy.Compress_management(workspace)
except:
    arcpy.GetMessages()


texts:
I am trying to download the polygons for the Umweltzone Heidelberg from OpenStreetMap using QuickOSM in QGIS. I queried by name as shown in the screenshot:

The corresponding OSM script looks like this:

<osm-script output="xml" timeout="25">
     <id-query {{geocodeArea:Germany}} into="area"/>
     <union>
          <query type="node">
               <has-kv k="name" v="Umweltzone Heidelberg"/>
               <area-query from="area"/>
          </query>
          <query type="way">
               <has-kv k="name" v="Umweltzone Heidelberg"/>
               <area-query from="area"/>
          </query>
          <query type="relation">
               <has-kv k="name" v="Umweltzone Heidelberg"/>
               <area-query from="area"/>
          </query>
     </union>
     <union>
          <item />
          <recurse type="down"/>
     </union>
     <print mode="body" />
</osm-script>

The query runs successful however no result is returned.

I also tried to download by relation id running this query:

rel(2676846); >; out meta;

This downloads two points - nothing more.

How can I download the polygons for the Umweltzone Heidelberg?


texts: Query layer in ArcGIS Desktop with PostgreSQL not querying data?

texts:
After some days I was finally able to import OSM data to a PostgreSQL database and connect it with ArcGIS 10.2.2

I can load every piece I want by simple clicking "Add data ..." however a SQL query layer always delivers the complete data rather than selecting the rows.

Here is what I set up first:

And that's what remains after ArcGIS finished loading the layer:

I don't understand why it keeps transforming the query expression.

In pgAdmin it shows just fine:

I am really clueless in this matter and ArcGIS is just so slow loading all the data over and over again (instead of QGIS, where the query is no problem).


texts: Setting band 4 of raster to alpha with GDAL?

texts:
I have a raster which has 4 bands, however band 4 is not being interpreted as an alpha channel. When I run gdalinfo I get the follow information about the bands:

Band 1 Block=360000x1 Type=Byte, ColorInterp=Gray
Band 2 Block=360000x1 Type=Byte, ColorInterp=Undefined
Band 3 Block=360000x1 Type=Byte, ColorInterp=Undefined
Band 4 Block=360000x1 Type=Byte, ColorInterp=Undefined

I would like it to look like this:

Band 1 Block=360000x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=360000x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=360000x1 Type=Byte, ColorInterp=Blue
 Mask Flags: PER_DATASET ALPHA
Band 4 Block=360000x1 Type=Byte, ColorInterp=Alpha

I am able to do this in QGIS by setting the band 4 to transparent and then rendering out the raster.

Is there a way of changing the "ColorInterp" in GDAL as I want to automate the process?


texts: Where does ArcMap GUI get font size and docking guide positions from?

texts:
In the screenshot below, ArcMap 10.1 appears to be making some arbitrary choices about the font to use for:

The Layers data frame name in the Table of Contents (larger font than expected)
The Main Menu pulldowns (smaller font than expected)

Can anyone tell me where to find a setting that will enable me to make these the same size as the font used for the Main Menu items, please?

I am using Windows 7 64bit and I know I messed around with Large Icons on the Options tab of Customize | Customize a while back.  I think the small font pulldowns started when I turned them back to normal.


texts: QGIS attribute table shows attributes of only one polygon at time

texts: I have multiple shapefiles each with many polygons. I am used to the attribute table listing them all in rows and columns. In this case, all attribute tables show me the attributes of a single polygon which I can select from a list on the left (see screenshot). How do I change my attribute table to the standard view? My QGIS version is 3.14


texts: Easytouse Android app to visualise offline maps

texts:
I need to propose an application android, ultrasimple of use which proposes to visualize various layers offline (french Scan25, Orthophoto). The idea would be to create mbtiles (with QGIS, Mobac, Mapcache, or equivalent) and prepare it for each employee, a file that corresponds to its intervention department, knowing that it must be 100% offline.

Do you have one (or more) applications to recommend that respond to that? The employees for some are novices in computer science, and it would have to be hypersimple to use (OruxMaps, Osmand, Qfield, etc.) offer too much.
Here, we would just like to offer the display of these funds, possibly choose the layers / transparencies, focus on the position of the GPS, possibly register a GPS point and that's it.

The better would be that we can offer a link to each employee who can recover the application and the file that contains the data (mbtiles format or equivalent), if possible free / not paid because there are many people to equip .

Would anyone have any advice?


texts: Add all vertexes X & Y value to attribute table

texts:
Is it possible to automatically extract every pair of (X,Y) coordinates of every vertices for a polygon and add it to its attribute table?

For example: if a polygon has 7 vertices I want 14 new fields in the attribute table: X1,Y1,X2,Y2,X3,Y3,X4,Y4.....
Or if it has 4 vertices I want the fields X1,Y1,X2,Y2,X3,Y3,X4,Y4 to contain their coordinates and the other ones (X5,Y5,X6,Y6....) to be empty.

Is there a tool / script that does such a thing?


texts: GPS Elephant Tracking II: Is ArcGIS basemap imagery orthorectified?

texts:
-- For my original question and pictures of elephants wearing GPS equipped anklets, go here. --

I'm experimenting with the use of GPS anklets to track the movement of zoo elephants.  This has proved endlessly complex, but I'm focused on one particular question now: is ArcGIS basemap imagery orthorectified?

My elephants wear GPS anklets that collect data points every 5 seconds.  My goal is to cleanup this raw data as best as I can, without going overboard.  Some of the steps I've taken are applying a variety of filters (NSAT>=4, HDOP<2, min and max movement thresholds based on elephant speeds) and keeping track of the times when the elephants are inside structures so that I can remove those points from the dataset.  When I map the remaining points, some are still falling outside of the exhibit boundaries.

Originally, I had planned to map the data on the basemap, draw a polygon on top of the elephant exhibit, and use the clip function to remove points falling outside the exhibit.  I've heard that I'm operating on the flawed assumption that the basemap provides an accurate view of the exhibit, when in fact, various distortions (camera tilt, local topography) may be skewing the imagery.  I've also read that ArcGIS buys it's imagery from various sources, all of which probably do some orthorectification, but that there is no way to know the accuracy of the map at any particular place.  On the other hand, my research may be operating on such a small scale that these distortions, if they exist, are unlikely to cause any problems.

Any solution would need to be easily implemented and widely applicable, as this data is being collected at zoos across North America.

Thanks,
Matty


texts: model builder expressions code block and definitions

texts: I am having difficulty wrapping my brain around the point of the expression vs the code block. What am I suppose to be defining? I copied the example from calculate value in the ArcGIS help website almost verbatim and yet it says my definition is wrong. More specifically it says a syntax error in line 1.


texts: Raster visualisation in equal area projections

texts:
This is a rather conceptual question that i haven't been able to address yet. Let's say I have a raster file projected in an equal area projection (say Mollweide), with a cell resolution of 100km x 100km. Since the projection preserves the area along the parallels, all cells should have the same area.

Now, as we know from the basic principles of cartography, no projection can preserve both area and shape of features when converting a sphere to a flat surface. This can be easily illustrated, if we draw a line of 100 km projected length at the equator on a map in equal area projection (but also on a map in orthographic projection), and a line of the same projected length towards the poles, these two lines will appear to have very different lengths in terms of screen pixels (check it out in google earth!).

Given this knowledge, if we try to visualise the above mentioned raster either in Python or in R using the appropriate packages (see scripts below), the cells have the same screen pixel dimensions regardless of the latitude. However, I would expect that, since each cell has the same dimensions (100km x 100km), cells towards the poles should be thinner and rather elongated (see also the attached snapshot from a recent scientific publication) in order to preserve area. So my question is, what is the appropriate way of visualising raster cells in equal area projections?

Code in Python:

import matplotlib.pyplot as plt
from osgeo import gdal
import cartopy
import numpy as np

fig = plt.figure(figsize=(15, 8))

proj = cartopy.crs.Mollweide()
ax = plt.axes(projection = proj)

raster = gdal.Open('test.tif') # a raster at with cell resolution 100x100km
band = raster.GetRasterBand(1)
nodata = band.GetNoDataValue()
array = band.ReadAsArray()

array[array == nodata] = np.nan
img = ax.imshow(array, origin='upper', cmap='jet')

ax.get_xaxis().set_visible(True)
ax.get_yaxis().set_visible(True)

plt.show()

Code in R:

library(raster)
ras = raster('test.tif')
plot(ras)

Example figure with stretched cells towards the north pole (source):


texts: Torque.js is not working: L.TorqueLayer is not a constructor

texts:
I've been trying to create a carto map with a torque layer as indicated in the Torque.js documentation, but I keep getting the same error:

Uncaught TypeError: L.TorqueLayer is not a constructor

<script>
  var CARTOCSS = [
    'Map {',
    '-torque-time-attribute: "date";',
    '-torque-aggregation-function: "count(cartodb_id)";',
    '-torque-frame-count: 760;',
    '-torque-animation-duration: 15;',
    '-torque-resolution: 2',
    '}',
    '#layer {',
    '  marker-width: 3;',
    '  marker-fill-opacity: 0.8;',
    '  marker-fill: #FEE391; ',
    '  comp-op: "lighten";',
    '}'
  ].join('\n');

  var torqueLayer = new L.TorqueLayer({
    user: 'your_username',
    table: 'your_table_name',
    cartocss: CARTOCSS
  });
  torqueLayer.addTo(map);
  torqueLayer.play()
</script>

How can I make L.TorqueLayer work?


texts: ESRI JS API - Unable to register Basemap Gallery (ID is already registered)

texts:
I am attempting to add a BasemapGallery to my web application. I prefer to create my layout dynamically using constructors such as:

layout.cpRight = new ContentPane({
    region:"right",
    ...
    ...
});

layout.borderContainer.addChild(layout.cpRight);

This is functioning fine and my layout is placed as it should.

The probleme is When I attempt to add basemapGallery with the bellow code:

maps.basemapGalleryParams = {
    showArcGISBasemaps:true,
    map: maps.map,
};

maps.basemapGallery = new BasemapGallery(maps.basemapGalleryParams,layout.cpRight);
maps.basemapGallery.startup();

I get the following Error Code: init.js.160

  Error: Tried to register widget with id==dijit_layout_ContentPane_0
  but that id is already  registered(…) "in domReady callback" "Error:
  Tried to register widget with  id==dijit_layout_ContentPane_0 but that
  id is already registered

I originally declared ids inside of the widget constructors parameter but removed when I got this error. The results where the same.

A few web searches indicate this can be a fairly common problem with some "easy" solutions but I have not found a solution where the layout is declared programaticly  as I would like to. It seems many defer to declare components in HTML.

Side note:
I break out applications into multiple scripts thus I create arrays to ensure certain variables are global across all scripts. Hence why in the scripts above, things are not being declared with the var keyword, but rather directly appending to the array.


texts:
I'm on a Windows 7 machine with QGIS 2.8 installed via OSGeo4W. I'm developing several processing scripts that make use of a few common pieces of code. Rather than maintain the code in each individual script file, I'd like to import the common functions into each script.

These scripts may be shared between multiple users so I'd like to avoid having code stored in multiple locations. My initial thought was to save _foo.py within the .qgis2\processing\scripts folder on my user account and then import via

import _foo

However, Processing doesn't set the scripts directory as the CWD of the python interpreter so the utility script doesn't get picked up.

Is there a way to achieve my goal without turning this into a full-blown QGIS plugin?


texts:
Based on my previous question, answered here Can not get the connecting line to align with label in desired place in QGIS 3. Now I have a new issue.
As can be seen in the image included, I have manged to get leaders from my points to my labels. The problem is, at the label they are anchoring to the bottom left hand corner of the label, which for this image just does not look right. My preference is to make the leader line anchor to the top left hand corner of the label text or the centre of the label, again at the top of the label.
I have tried to alter the horizontal (to Centre) and vertical (to Half) as suggested in a similar question posed by someone else but this has not worked.
How can I change the anchor point to either of the above options?


texts: TIF not showing UTM coordinates in world file

texts:
I have a lot of tif files (ortho imagery that came with Lidar data). Each tile has a .tif and .tfw but the world files (.tfw) have line 5 and 6 as numbers that aren't what I would expect as I was expecting the files to be in UTM ZONE 10.

Is there anything I am missing that usually comes with Ortho Imagery / Lidar?
Or do I need to contact the supplier and ask them to re-export as a geo-referenced file.
Or maybe this is just coordinates on a local site datum and that is it.

Below is an example of on the of the world files. The list number is the line number.

0.200078983
0
0
-0.200078983
-1813.10469
4172.72779

I guess the question is can you have ortho-rectified imagery on a local site datum? I thought maybe it would have to be real-world coordinates?


texts:
I am using GeoExt, GeoServer, and PostgreSQL 8.2. I have a shapefile of trees and want to show the total number of trees to the user...

Is there any feature in GeoExt like feature count which can make my work easy such as a distance and area calculation actions in a tool bar?


texts:
I am new to WebODM and I am trying to generate DSM and Orthomosaic images using default settings in WebODM.

The Orthomosaic image is correct by providing GCP but the DSM file is having inaccurate elevations when compared to google earth.

I am not sure what went wrong or is there any way/configurations in WebODM to have a better DSM generated.

I don't know what parameters should I use for better results so I keep the input parameters Default or DSM + DTM.
I'm attaching a GCP table, screenshot of DSM, Orthophoto, and input parameters, so can you please suggest to me what changes should I made to input parameters for getting better output.

GCP table:
+proj=utm +zone=44 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
357783.55   3120786.54  0   5178.90 1405.57 DJI_0001.JPG    gcp16
357783.55   3120786.54  0   5162.22 2005.07 DJI_0002.JPG    gcp16
357783.55   3120786.54  0   5142.58 2778.84 DJI_0003.JPG    gcp16
357804.03   3120770.86  0   4262.03 1258.89 DJI_0001.JPG
357762.67   3120667.30  0   4948.59 2184.14 DJI_0007.JPG
357777.35   3120672.63  0   5355.83 1863.45 DJI_0007.JPG
357754.61   3120678.37  0   4503.34 1714.50 DJI_0007.JPG
357760.16   3120490.52  0   3624.37 1034.95 DJI_0016.JPG
357731.51   3120529.78  0   4986.50 1891.50 DJI_0016.JPG    gcp15
357743.91   3120519.01  0   4411.38 1939.29 DJI_0016.JPG


texts:
How to put the field from arcpy.GetParameterAsText() into the expression here:
CalculateField_management(in_table, field, expression, {expression_type}, {code_block})

Now it is like "!Name!.replace(\"о\",\"r\")"
but here the field !Name! is already firmly defined and isn't depend on arcpy.GetParameterAsText()


texts: Using AndNav tiles storage as input for OpenLayers.TMS

texts:
Is there a possibility to use map tiles from AndNav (created wit Mobile Atlas Creator) that are in format:

Mapnik\12\2241\1369.png.andnav (I know that 12 is the zoom level but don't have idea what is 2241 and 1369...)

as an input for OpenLayers.TMS layer?

I found information that I could use local tiles with TMS: http://trac.osgeo.org/openlayers/wiki/UsingCustomTiles

Unfortunatelly, I can find any place where AndNav map format is documented and reading sources didn't helped me too much.

Would anyone could help me?


texts: QGIS cluster character symbol visibility based on cluster size

texts:
I want to display the cluster size with a condition : I don't want to display this information for small clusters.
The expression is working fine but I don't know I can add else...

case
when @cluster_size > 30 then @cluster_size
end

Is my case, cluster size of the clusters with a size > 30 are displayed. BUT cluster <30 are labelled with a A letter AND I want to display NOTHING.

Not working :

case
when @cluster_size > 30
then @cluster_size
else @cluster_size = NULL
end

Returns me 0 but not empty :

case
when @cluster_size > 30
then @cluster_size
else @cluster_size = ' '
end

If I remove A letter from the assistant, nothing is displayed.


texts: Why is normalize_height() function aborting in lidR?

texts:
I'm trying to run normalize_height() on 328 las files. Maybe or maybe not relevant: I'm not running using LAScatalog because I couldn't get that to run at all. Anyhow, I'm running using a parallelized foreach loop. Seems to work fine. 303 of the tiles were properly normalized and written to disk. For others I get this output (looping/parallelization not relevant):
library(lidR)
las <- readLAS("1715.las")
las <- normalize_height(las, tin())
#> Errors running normalize_height: "270 points not normalizable. Process aborted.

When I run las_check() on this file the suspicious results are:

checking gpstime incoherences: 405091 pulses (points with the same gpstime) have points with identical ReturnNumber
checking negative outliers: 177777 points below 0

My question is: how do I track down this problem? All I could think of was removing duplicates, which I did with no success. Do I need to remove the degenerated ground points too? If so, how to do that? The point cloud that triggered the error can be downloaded here and looks like:


texts: Measure of "Matchingness"/Correlation between polylines

texts:
I am working with two layers of polylines that represent different modeled/representative approaches to summarizing aircraft radar tracks.

I am curious if, beyond just visually inspecting them, if there is an appropriate quantitative way of looking at how similar or "related" they are in terms of their direction and turns, etc, rather than something basic like looking at the amount of parts along the polyline layers that do and do not intersect.

I want to say there is a way at looking at a correlation between the two, but something about that seems rather complicated (which I am willing to read more about).  I also am not sure which GIS/spatial data software would be the most useful for this sort of task (i.e. QGIS vs GeoDA vs ArcMap).

Therefore, I am wondering two fold:
1) What sort of statistical test/method would be most appropriate?
2) Which spatial platforms would be best for carrying it out- does someone know specifically which packages/features/commands would be most useful within each platform?


texts: Masking/removing polygons from shapefile

texts:
I have a shapefile containing the road traffic noise of a city.

The shapefile has 5 columns, where 'noise_in' (noise interval) ranging from 1-5 is included.

I want to remove/mask the road noise from a specific area.

My problem is that the noise intervals covers more than just this area, e.g. if I remove a polygon from the file, covering the area of interest, it will also remove the noise from an adjacent area, which I'm not interested in.

The attribute table looks like this:

The road noise for the area looks like this:

I would like to remove all the road noise from this area, meaning all 5 noise intervals. Unfortunately, I have no available date for recalculating the road noise for the area.

I don't have a polygon for each area.

Here's an example of "Identify Features":

The red highlighted areas are two different polygons covering large areas of road noise, and I'm only interested in the big blue "vein" in the middle.

I'm using QGIS 1.8 Lisboa for MAC

Thanks in advance