texts: ValueError when reading a raster band as a numpy array

texts:
Here's the problem shown in Python shell:

In [32]: from osgeo import gdal

In [33]: geotiff = '/Users/calvin/work/luc/proj_public/media/GLO
/Users/calvin/work/luc/proj_public/media/GLOBCOVER_L4_200901_200912_V2.3.color.tif
/Users/calvin/work/luc/proj_public/media/GLOBCOVER_L4_200901_200912_V2.3.color_orig.tif
/Users/calvin/work/luc/proj_public/media/GLOBCOVER_L4_200901_200912_V2.3.tif
/Users/calvin/work/luc/proj_public/media/GLOBCOVER_L4_200901_200912_V2.3_CLA_QL.tif

In [33]: geotiff = '/Users/calvin/work/luc/proj_public/media/GLOBCOVER_L4_200901_200912_V2.3.color.tif'

In [34]: src = gdal.Open(geotiff)

In [35]: band = src.GetRasterBand(1)

In [36]: band
Out[36]: <osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x10632b8f0> >

In [37]: ar = band.ReadAsArray()
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/Users/calvin/work/luc/<ipython-input-37-d7dde170464f> in <module>()
----> 1 ar = band.ReadAsArray()

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/osgeo/gdal.pyc in ReadAsArray(self, xoff, yoff, win_xsize, win_ysize, buf_xsize, buf_ysize, buf_obj)
    893         return gdalnumeric.BandReadAsArray( self, xoff, yoff,
    894                                             win_xsize, win_ysize,
--> 895                                             buf_xsize, buf_ysize, buf_obj )
    896
    897     def WriteArray(self, array, xoff=0, yoff=0):

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/osgeo/gdal_array.pyc in BandReadAsArray(band, xoff, yoff, win_xsize, win_ysize, buf_xsize, buf_ysize, buf_obj)
    225                                     buf_xsize, buf_ysize, datatype )
    226         ar = numpy.fromstring(band_str,dtype=typecode)
--> 227         ar = numpy.reshape(ar, [buf_ysize,buf_xsize])
    228
    229         return ar

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in reshape(a, newshape, order)
    169     except AttributeError:
    170         return _wrapit(a, 'reshape', newshape, order=order)
--> 171     return reshape(newshape, order=order)
    172
    173

ValueError: total size of new array must be unchanged


texts: How much of these USGS topographical contour maps are artistic license?

texts:
I was looking at the below map and got to wondering how much of all of those minute squiggling contour lines are real detail and how many are just artistic flair:

I could imagine that satellite laser/radar imaging could generate this image computationally today, but this map was printed in 1956, back when mapmakers were still etching plates by hand.

Did people actually scour the entire United States measuring elevation at this insanely fine level of precision to produce these maps? To what extent are the contour lines inferred from coarse-grained (maybe aerial) observations like "there's a ridge over there, and it's kind of jagged, so we'll make it squiggle a few times" versus "there are exactly nine squiggles between these two creeks 100 meters apart?"


texts: Insert point layer as input for GRASS viewshed in QGIS 3.2 modeler

texts:
I’m building a QGIS (3.2 Bonn) tool for my colleagues using graphical modeler.
The tool should get 3 inputs: DEM raster layer, vector layer of points, vector layer of a point of interest (POI). The POI should be different each time the user executes the model. One of the processes in the model should be to calculate viewshed using the DEM raster and the POI (vector layer). I saw that the only tool available in QGIS 3 for modeler is GRASS r.viewshed which gets as input raster DEM and set of xy coordinates to identify the viewing position. Is there a way to insert instead of the xy coordinate the POI, meaning a point layer? or somehow calculate automatically (in the modeler) the xy coordinates of the POI the user choose and then extract the xy information from the layer and insert it to the viewshed tool?

My workaround so far is a bit cumbersome, in the modeler I calculate xy coordinates using add_geometry_attributes and use the field calculator to concatenate the x,y and insert it to a new field. Then, I manually copy the new field value (e.g., 115.45, -22.455  [EPSG:4326]) and manually paste it in the viewshed tool.

I saw some answers for similar question, but they are relevant to older versions of QGIS, for example this but it uses table parameter that does not exist in QGIS 3 modeler, I also tried this but I can't see the script in the algorithm section in the modeler.

I wrote some basic QGIS tools in the python editor in QGIS 2 but it seems that in the new version things have changed and the is a lack of documentation.


texts:
I recently downloaded and installed  GeoServer 2.6.2.
One of my other colleagues asked me to add an SQL store for a proof of concept we are doing.

When navigating to the "add new store" page I could not find the SQL store.
I have followed the instructions in the documentation, but cannot seem to get it working.

I have tried GeoServer does not show sql data store?, but I still cannot seem to get it configured.

I am new to web mapping.

My system:
 - Win7 Pro.
 - jre-8u40-windows-i586
 - GeoServer 2.6


texts: QGIS unable to join attribute by location - attribute table missing values

texts:
I have a shapefile, clean_intersection, representing a geographic area. I'm trying to add additional columns to the attribute table showing the evapotranspiration (ET) for each element of my map, clean_intersection.

I used the "simplify geometry" tool to simplify clean_intersection enough to handle in Excel, where I created a .csv file to add a couple columns to the attributes table. Then I used "add delimited text layer" to return the .csv file to QGIS now as "clean_simple_data." Clean_simple_data is a shapefile with the data I want but with simplified geometry polygons. I'm now trying to join clean_simple_data with clean_intersection to generate a shapefile with detailed boundaries (like clean_intersection) but with the additional columns from the clean_simple_data attribute table.

To accomplish this join, I try:

This creates a joined layer. However, the attribute table of the joined layer has many repeated values:

The "ICA" value of 6.633 is correctly paired with the ET value of 1.12, but it is assigned to too many rows. Here's a screenshot of the Excel .csv file showing how there should be many different values instead of a few repeated values in the ICA and ET columns.

Originally I had many more columns of data that I was attempting to join but I pared the .csv file down which seemed to help but did not resolve the issue. I have also tried many different variations of the geometric predicates available in the "Join attributes by location" dialogue, to no avail. I am running QGIS 2.18.

Does anyone have any insight?


texts: Calculating distance from line to multiple points in ArcGIS Desktop?

texts:
I have road maps (line features) and point features representing forest plots. I created buffer around the roads. Now I want to know the distance from line(road) to the points within this buffer.

How can I do that?


texts: Census.gov usually makes available the shapefiles containing school zone boundaries. Due to the US government shutdown, census.gov is down. Does anyone know where to download a recently archived copy?


texts:
this is a rather simple and trivial question, but before I contact my supervisor I want to get sure. I have a .tif image that says the resolution is 4.5 x 4.5 mm. The output of gdalinfo is: Pixel Size = (0.004545650000000,-0.004545649999975). When I look at the image it looks very pixelated.

In my view a 4.5 mm resolution means that one pixel displays 4.5 mm x4.5 mm of the real view?

EDIT

Size is 53243, 11916
    PROJCS["WGS 84 / UTM zone 32N",
        GEOGCS["WGS 84",
            DATUM["WGS_1984",
                SPHEROID["WGS 84",6378137,298.257223563,
                    AUTHORITY["EPSG","7030"]],
                TOWGS84[0,0,0,0,0,0,0],
                AUTHORITY["EPSG","6326"]],
            PRIMEM["Greenwich",0],
            UNIT["degree",0.0174532925199433],
            AUTHORITY["EPSG","4326"]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["latitude_of_origin",0],
        PARAMETER["central_meridian",9],
        PARAMETER["scale_factor",0.9996],
        PARAMETER["false_easting",500000],
        PARAMETER["false_northing",0],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AUTHORITY["EPSG","32632"]]
    Origin = (357377.652999999991152,5610076.298999999649823)
    Pixel Size = (0.004545650000000,-0.004545649999975)
    Metadata:
      AREA_OR_POINT=Area
    Image Structure Metadata:
      INTERLEAVE=PIXEL
    Band 1 Block=256x256 Type=Byte, ColorInterp=Red
      Unit Type: metre
    Band 2 Block=256x256 Type=Byte, ColorInterp=Green
      Unit Type: metre
    Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
      Unit Type: metre

    Corner Coordinates:
    Upper Left  (  357377.653, 5610076.299) (  6d59' 0.91"E, 50d37'30.36"N)
    Lower Left  (  357377.653, 5610022.133) (  6d59' 0.99"E, 50d37'28.61"N)
    Upper Right (  357619.677, 5610076.299) (  6d59'13.22"E, 50d37'30.57"N)
    Lower Right (  357619.677, 5610022.133) (  6d59'13.30"E, 50d37'28.82"N)
    Center      (  357498.665, 5610049.216) (  6d59' 7.11"E, 50d37'29.59"N

)


texts: Pixel depth and raster output won't align with the original

texts:
The problem

Need to convert a raster .tif tile (DOWNLOAD) from WGS84 to GDA94 Albers;
Need to convert its stored pixel info from Float32 to 8 bit UNSIGNED;
Need Compression set to "NONE";
Need to multiply it by 40; and
Need to set its pixel size to 30x30 metres.

My steps

1 - Project the raster tile to GDA94 Albers, stating "30" as Output cell size: ArcMap -> Project raster tool;

2 -  Recalculated the pixel values multiplying by 40: QGIS -> Raster calculator (Syntax = "tile_name" * 40);

3 - Convert the pixel depth from Float32 to 8-Unsigned: ArcMap -> copy raster tool.

My results

As you can see from the picture below, I have sampled some points and retrieved the info from the original tile and compared to my output raster:

After multiplying by 40, the result I SHOULD expect should be equal to the "OrigX40", however, my results are shown in the "New" column and there's a substantial difference.

In addition, is it normal that after reprojecting the final output it shifts the cells' orientation from tilted (pic above) to non-tilted (pic below)?

This shifted array of cells also contributes to the fact that when I sample a point that is on the fringe between 2 cells of the final output, obviously it does not match with the expected Float value (divided by 40) of the original tile. Is this normal or should the final output maintain a shifted array-like pixel grid, same as the original?

My tools

I have used a combination of QGIS 3.6 Noosa and ArcMap 10 (no Spatial Analyst extension -> no raster calculator!).

My expectations

My expectations are that the final output should align with the original layer and every cell stores the same info the original cells store but times 40.

What am I doing wrong?


texts:
I have a large list of raster images (>600) I would like to automatically georeference to a source image.

I can do this in ArcGIS 10.2 using the Auto Registration tool for each image but looking for a way to iterate it in ModelBuilder/python script.

Is there an arcpy command or similar for the auto registration tool?


texts: GRASS: Can v.buffer use geodesic distance like r.buffer does?

texts:
Question: In GRASS GIS, can v.buffer use geodesic distance like r.buffer does? If v.buffer can't use geodesic distance, is there a way to do vector geodesic buffering (ideally in GRASS, but I can switch to something else if necessary)

I'm a GIS newbie, and find this really frustrating:

The r.buffer man page clearly states that geodesic distances are used for latitude/longitude coordinates:

  r.buffer measures distances from center of cell to center of cell
  using Euclidean distance measure for planimetric locations (like UTM)
  and using ellipsoidal geodesic distance measure for latitude/longitude
  locations.

The v.buffer man page doesn't say this, but it doesn't say the opposite either. In fact, the word "Euclid" appears nowhere in this document.

I've google'd around (maybe using the wrong phrases, and the dot in "v.buffer" and "r.buffer" may confuse search engines), and vaguely remember seeing something that says GRASS GIS v.buffer can not create geodesic buffers regardless of what options you use. However, I'm not sure I saw this and didn't bookmark the page since I was looking for something else at the time.

I tried using "v.extrude" and then "v.buffer" on the resulting geometry, but that didn't work either.

My "goal" is to solve https://earthscience.stackexchange.com/questions/14656/how-to-calculate-boundary-around-all-land-on-earth just for "fun". I'm aware of https://oceancolor.gsfc.nasa.gov/docs/distfromcoast/ but want to see if I can either re-create it or do a "better" job by using a higher resolution coastline from https://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/gshhg-shp-2.3.7.zip (ie, using the "f" level resolution instead of the "c" level that the previous URL uses).


texts:
I’m working on some maps for the region, and assessing available datasets. So far I have found a few sites with interesting data, but not many. Does anyone have any suggestions about where I could look? Mostly looking for vector data.

Thanks!

Edit: I am focusing mostly on human data, such as religion, health, etc, which could include some physical features, such as buildings. Overall though, I'm just trying to see what is out there.


texts: Convert a GeoJSON polygon to a Cesium.Plane object

texts:
I have a polygon (planar) in a GeoJSON feature.
I would like to convert it to a true Cesium.Plane object.
How could I achieve that?
In my current code, the feature is only a dataSource object:
const viewer = new Cesium.Viewer("cesiumContainer");

const GeoJSON = {
    "type":"FeatureCollection",
    "name":"MyPolygon",
    "crs":{"type":"name","properties":{"name":"urn:ogc:def:crs:OGC:1.3:CRS84"}},
    "features":[{"type":"Feature","properties":{},"geometry":{"type":"Polygon","coordinates":
      [[[10.746500,48.31470,500],
        [10.747500,48.31570,500],
        [10.747038,48.31590,550],
        [10.746038,48.31490,550],
        [10.746500,48.31470,500]]]}}]
};

const loadData = async () => {
    try {
        const dataSource = await Cesium.GeoJsonDataSource.load(GeoJSON, {
                stroke: Cesium.Color.BLUE,
                fill: new Cesium.Color(0.2, 0.2, 1, 0.2),
                strokeWidth: 3
        });
        viewer.dataSources.add(dataSource);
        viewer.zoomTo(dataSource);
    }
    catch (err) {
        console.log("Error: ", err);
    }
};
loadData();

https://sandcastle.cesium.com/


texts: Change cellsize of ArcInfo Grid ASCII data

texts:
I have a ArcInfo Grid ASCII file of Turkey with cellsize 0.000833333333 (3 arc second ~90m for equator and ~83m for Turkey, I've learned so far), Can I just simply change the header value to 83m or is there any proper way to change cell size to an integer number. How can I do that?

ncols        24001
nrows        12001
xllcorner    24.999583987092
yllcorner    34.999583575447
cellsize     0.000833333333
NODATA_value -9999


texts: QGIS 3.14: new attribute sorting feature

texts:
The 3.14 Changelog includes this new feature:
"Allow sorting attribute table by field not listed in the table"
The description states: "Instead of using the same data model for the displayed and the sorting columns, two data models are now used. They use the same API / base class as they are very similar."
This new feature sounds intriguing, but I don't understand what the changelog is saying, nor can I find any other help on the web.
When and how would I use this feature?


texts: Georeferencing AutoCad Drawing in ArcGIS Pro

texts: I requested an AutoCAD file of a small town's streets/buildings to try out georeferencing the drawing as practice. Needless to say, this drawing does not have a coordinate system specified, but the measurements in the drawing are pretty accurate and are in FT. However, when I specify the coordinate system in AutoCAD it is telling me that the measurements are the same, but on the basemap it covers a third of the state. I have the same thing happen in ArcGIS Pro except it gives me the actual measurements (around 30x larger than normal).  I was hoping there would be a way to easily scale it down to size, but that doesnt seem to work. I've tried using the georeference tool in ArcGIS Pro and that does help, but it still ends up being about 30% bigger than the basemap regardless.


texts: Joining point layer with line layer using coordinates

texts:
I am using QGIS 3.14, and I want to represent a road network.

I had a polyline layer (3257 entities) as a background. For each one I have calculated x_start, y_start, x_end, y_end coordinates, length of each line and other elements.

I have also created a point layer with every start point (3257 entities) in which I also calculated x, y for each point. I did the same for the 3257 end points.

These point layers have also altitude information which I extracted from a DEM. I want to transfer this altitude information to line layer.

I tried to use JOIN ATTRIBUTES BY FIELD VALUE and used x_start (line layer) and x (start point layer) as "key", and did the same with x_end (line layer) and x(end point layer), in order each line afterwards to have the correct altitude of its startpoint and endpoint. But result is correct for some only of the lines (I checked it via excel).

What have I done wrong?


texts: Why is ESRI.ArcGIS.Desktop.AddIns.ComboBox configured with a number of W's?

texts:
I'm using a ESRI.ArcGIS.Desktop.AddIns.ComboBox. I configure the size of the combobox in config.esriaddinx with the following xml:

    <ComboBox id="my_id"
              class="MyComboBox"
              message="Everything is awesome!"
              caption="Everything is awesome!"
              tip="Everything is awesome!"
              category="Add-In Controls"
              sizeString =     "WWWWWWWWWWWWW"
              itemSizeString = "WWWWWWWWWWWWW" />

It work fine, but it really baffles me that the width is configured with a bunch of 'W's.

Can anyone explain the reasoning behind it?


texts: Clustering points into *n* number of groups of equal counts with ArcGIS

texts:
I have to  divide points into groups. These groups will be used by agents to talk to property owners. The goal is to make the agent's job easy by grouping parcels that are near each other together, as well as divide the parcels into (almost) equal numbers so that the work is distributed evenly. The number of agents is variable: from 2 to 8.

I've explored various ArcGIS tools (Grouping analysis and Distriction), but none seem to suit my need.

Is there a reliable method to do this task?


texts: Tilestache Multiple Grids results in PostGIS error

texts:
I've got a PostGIS database which Mapnik and Tilestache pull from. The mapnik xml has multiple layers as I'm pulling from different tables depending on the zoom - and these different tables have different column names.

The issue is that I'm trying to get around is that it seems using multiple grids within tilestache requires the column names to be the same.

Here's the tilestache cfg:

"grid":
        {
            "provider":
            {
                "name": "mapnik grid",
                "mapfile": "/usr/local/share/maps/style/OSMBright/OSMBright.xml",
                "layers":
                [
                    [0, ["name", "employees"]],
                    [1, ["name", "employees"]],
                    [2, ["id", "employees"]],
                    [3, ["company_name", "employees"]]
                ]
            }
        }

Zooms at layers 0 and 1 work because they have similar column headers. It fails starting at layer 2:

  RuntimeError: Postgis Plugin: ERROR:  column "name" does not exist
  LINE 1: ...LECT ST_AsBinary("geom") AS geom,"employees","id","name"
  FRO...
                                                               ^ in executeQuery Full sql was: 'SELECT ST_AsBinary("geom") AS
  geom,"employees","id","name" FROM ( SELECT geom, id, employees
                FROM blkgrp_2014
              ) AS data
               WHERE "geom" && ST_SetSRID('BOX3D(-10693846.00520928 5028944.964938281,-10684062.06558878 5038728.904558787)'::box3d, 3857)'

For some reason it's looking for the column "name" even though my config has specified to look for "id" - so tilestache is obviously progressively combining the search.

Is there an obvious solution to this? Or do I just add a null value column to each table so the column names match? Seems hacky...


texts:
I have created many .gpx datas with my Garmin device, but due to the amount of data, I have imported the mapping data from my device several times (also because I didn't want it to be lost)...

As a result I have now many double or even 6 times the same point (including tracks, which I displayed as points also, so that I have the elevation) after I saved all the gpx-files as shapefiles and then merged them.

Is there a way to automatically erase the duplicates with QGIS, or if not with any other free software working on Ubuntu 12.4?

I very much appreciate your help!


texts:
I have a shapefile (sf) object made up of polygons with polygon characteristics stored in a dataframe. This object is called tiles. Within the dataframe there are 12 variables indicating number of pixels that transitioned to zero forest cover by year called loss2001, loss 2002, ... loss2012. Each variable is indexed by a unique identifier ID.
library(sf)
tiles = st_read("ShapefileForestlossGroup1")

I want to reshape this shapefile dataframe so that I can produce annual spatial weights. The problem is that many reshaping options are best executed with a data table, but if I transform tiles into a data table, it is no longer an sf object.
As an sf object with a dataframe, I tried executing the following reshape function:
tiles2 = reshape(tiles,
        direction = "long",
        varying = list(names(tiles)[1:12]),
        v.names = "numPixelsDeforested",
        idvar = "ID",
        timevar = "Year",
        times = 2001:2012)

But R returns the following error:
Error in `.rowNamesDF<-`(x, value = value) : invalid 'row.names' length In addition: Warning message: In `[<-.data.frame`(`*tmp*`, , v.names, value = list(loss2001 = c(0L,  : provided 2 variables to replace 1 variables
Could these errors be related to the fact that this is an sf object dataframe, or is the syntax just wrong?


texts:
I keep trying to load Google Maps to start a project and I keep getting this error:

NotImplementedError: QgsPluginLayer.setTransformContext() is abstract and must be overridden

How do I fix it?


texts: Unexpected results when using ST_Intersection for polygon features

texts:
I've to create a spatial intersection of two tables (poly_a and poly_b) containing polygon features. Is there any way to remove unexpected linestring results?

CREATE VIEW poly_intersection AS SELECT
row_number() over() AS gid,
subquery.*
FROM
    (SELECT
    (ST_Dump(ST_Intersection(poly_a.geom, poly_b.geom))).geom::geometry(polygon, 4326)
    FROM poly_a, poly_b
    WHERE ST_Intersects(poly_a.geom, poly_b.geom) AND NOT ST_Touches(poly_a.geom, poly_b.geom)
    ) AS subquery;

ERROR:  Geometry type (LineString) does not match column type (Polygon)

I'm wondering about the line geometries because of using AND NOT ST_Touches(poly_a.geom, poly_b.geom) in the WHERE clause.


texts: Calculating distance between two sets of Decimal Degree coordinates?

texts: I have coordinates in decimal degrees e.g: point A:Lat,Lng(5.317101667,-2.023303333) and point B:lat,lng(5.317093333,-2.023318333). The altitude in Meters at Point A is 11m and the altitude and Point B is 33m, How can i accurately calculate the distance between these points? Note I know about Haversine, howeve,r these points are close to each other and I want an accurate way to calculate the distance between these two, would finding the distance between the points with d=sqrt(x2-x1)^2+(y2-y1)^2+(Z2-Z1)^2, be accurate? if so what units would the answer be in? I am looking to get the answer in metres.


texts: Re-projecting from Polar projection to EPSG:4326?

texts:
I simply want to project a polar stereographic raster (EPSG:3411) to Mercator (EPSG:3857) or WGS84 (EPSG:4326). The data I is from NOAA.

$gdalwarp -s_srs EPSG:3411 -t_srs EPSG:3857 northpsg.20141027 output.tif

The resulting raster is croped near the height of 75°N - which is exactly the opposite behavior of projecting into polar stereographic that zetah has

On the other hand the projection into WGS84 results in a output with a height of 0.

Do you have any suggestions how this could work?

I already contacted the GDAL community but didn't get any responses yet.


texts:
I was wondering if there is something like a wishlist for QGIS plugins?
The link on this page does not work anymore.
It would be quite useful to know what people need.


texts: How can I draw an outline around my raster classes

texts: I have a raster created from polygon shapes from a social survey of fishing areas in ArcMap. I have classified my raster to 3 classes High,Low,Medium. I want to draw a border or an outline around the raster classes so it is visible like a hollowed effect in three colours so areas that are High, Low, Medium are visible. How can I do this?


texts: How to find a shapefile's correct EPSG code

texts:
My company uses inhouse software in order to survey and we need to manually specify an EPSG code in order to set up a project and then add shapefiles. However, there is some confusion on my part because this needs to be an EPSG code and the numeric code for the projections I often cannot find as an EPSG code
As an example, today I found out that I messed up a project because I could not find an EPSG code for the NAD_1983_StatePlane_Maryland_FIPS_1900_Feet projection. I cannot use the ESRI or SR-ORG codes, so what would I do next? Reproject all the shapefiles into a similar coordinate system that does have an EPSG code, how do I choose which to use?


texts:
I have an INCIDENTS point layer. Staff use the point layer in ArcGIS Desktop for day-to-day visual analytical purposes, such as:

Panning around maps that include the layer
Visually comparing the INCIDENTS point layer to other GIS layers
General exploration of the data

Some of the points overlap or nearly overlap. I'm concerned that the amount of points in a given area is not being effectively communicated to users.

How can I better display how many points are in a given area?


texts: R: Voronoi tesselation from long/lat data?

texts:
I want to create a SpatialPointDataFrame and Voronoi Tesselation from my .csv data downloaded from EPA, available here: http://uloz.to/xboQc7cu/no2-calif-csv

I convert my data in SPDF format in 2 ways:

my.tab<-read.csv("no2_calif.csv", stringsAsFactors = FALSE, sep = "")
coordinates(my.tab) <- ~ Longitude + Latitude

or

my.tab<-read.csv("no2_calif.csv", stringsAsFactors = FALSE, sep = "")
my.tab.spdf<- SpatialPointsDataFrame(cbind(my.tab$Longitude, my.tab$Latitude), my.tab)

and use Voronoi tesselation published here: http://carsonfarmer.com/2009/09/voronoi-polygons-with-r/

voronoipolygons = function(layer) {
    require(deldir)
    crds = layer@coords
    z = deldir(crds[,1], crds[,2])
    w = tile.list(z)
    polys = vector(mode='list', length=length(w))
    require(sp)
    for (i in seq(along=polys)) {
        pcrds = cbind(w[[i]]$x, w[[i]]$y)
        pcrds = rbind(pcrds, pcrds[1,])
        polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
    }
    SP = SpatialPolygons(polys)
    voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
        y=crds[,2], row.names=sapply(slot(SP, 'polygons'),
        function(x) slot(x, 'ID'))))
}

Using either

my.tab.voro <-voronoipolygons(my.tab)

or

my.tab.spdf.voro <-voronoipolygons(my.tab.spdf)

I obtain back the error:

 Error in data.frame(x = crds[, 1], y = crds[, 2], row.names = sapply(slot(SP,  :
  row names supplied are of the wrong length

Please, where the problem could be? do I need to have my coordinates in XY format? if yes, please, how? just reproject my data?


texts:
Im running arcmap 9.1. i have the code that lists feature classes in a list box from the SDE. The users select from one list box and adds it to another list box. How do i get it to show the layer name not feature class name?? How can i make sure it keeps the same symbology as it has in the SDE? When i run this code it adds it as new layer with new symbology. I do not have the password for the SDE so it is extracting the feature classes by using the parent workspace of a selected layer in the Layers list.

Public Sub AddSDECoverageLayer()

Dim pFWorkspace As IFeatureWorkspace
Dim pFeatureWorkspace As IFeatureWorkspace

Dim pFLayer As IFeatureLayer

Dim pMxDoc As IMxDocument
Dim pMap As IMap

Dim pDataSet As IDataset
Dim pFeatureLayer As IFeatureLayer
Dim pFeatureClassName As IFeatureClassName
Dim pUnknown As IUnknown
Dim pWorkSpace As IWorkspace
Dim pParentWorkSpaceName As IWorkspaceName

 ' Get the selected layer
 Set pMxDoc = ThisDocument

 Set pUnknown = pMxDoc.SelectedLayer
  If Not TypeOf pUnknown Is IFeatureLayer Then
   Exit Sub
 End If
 Set pFeatureLayer = pUnknown

 ' Make sure the layer is from SDE and get the parent workspace
   Set pDataSet = pFeatureLayer.FeatureClass
   Set pFeatureClassName = pDataSet.FullName
   Set pFWorkspace = pDataSet.Workspace
  ' If Not pWorkSpace.Type = esriRemoteDatabaseWorkspace Then
  '  MsgBox "Not an SDE feature class"
  '  Exit Sub
  'End If
  'Set pDataSet = pWorkSpace
  'Set pParentWorkSpaceName = pDataSet.FullName

  Dim intX As Integer

 'Add each FeatureLayer to the map document
 'Create New Layer
 ‘ How do I keep the symbology the same as whats on the SDE?

 For intX = 0 To ListBox2.ListCount - 1

     Set pFLayer = New FeatureLayer
        Set pFLayer.FeatureClass = pFWorkspace.OpenFeatureClass(ListBox2.List(intX))

        pFLayer.Name = pFLayer.FeatureClass.AliasName
        Set pMxDoc = ThisDocument
        Set pMap = pMxDoc.FocusMap

        pMap.AddLayer pFLayer
Next

pMxDoc.ActiveView.Refresh

End Sub

The code below gets the feature class name and adds them to listbox.

Private Sub CommandButton4_Click()

Dim pDoc As IMxDocument
Set pDoc = ThisDocument

Dim pUnknown As IUnknown

Dim pFeatureLayer As IFeatureLayer
Dim pMxDoc As IMxDocument
Dim pSdeWorkspaceFactory As IWorkspaceFactory

Dim pWorkSpace As IWorkspace
Dim pParentWorkSpaceName As IWorkspaceName
Dim pDataSet As IDataset
Dim pFeatureClassName As IFeatureClassName

Dim pFWorkspace As IFeatureWorkspace

 'Get the selected layer
  Set pMxDoc = ThisDocument
  Set pUnknown = pMxDoc.SelectedLayer
    If Not TypeOf pUnknown Is IFeatureLayer Then
    Exit Sub
   End If
   Set pFeatureLayer = pUnknown

  ' Make sure the layer is from SDE and get the parent workspace
  Set pDataSet = pFeatureLayer.FeatureClass
  Set pFeatureClassName = pDataSet.FullName
 'Set pFWorkspace = pDataSet.Workspace
  Set pWorkSpace = pFWorkspace
    If Not pWorkSpace.Type = esriRemoteDatabaseWorkspace Then
    MsgBox "Not an SDE feature class"
  Exit Sub
  End If
Set pDataSet = pWorkSpace
Set pParentWorkSpaceName = pDataSet.FullName

Dim pEnumDSName As IEnumDatasetName
Dim pSdeDSName As IDatasetName

Set pEnumDSName = pWorkSpace.DatasetNames(esriDTFeatureClass)
Set pSdeDSName = pEnumDSName.Next
Set pFWorkspace = pWorkSpace

how do I get it to show layer name and not feature class name?

   While Not pSdeDSName Is Nothing
    ListBox1.AddItem pSdeDSName.Name
  Set pSdeDSName = pEnumDSName.Next

Wend

End Sub


texts:
In ArcMap; during the Import Symbology function you had the option to import the "Complete Symbology Definition". However, I can not find the option in Pro or any Python script that could solve this problem.

I've read that using the "Apply Symbology From Layer" geoprocessing tool would accomplish the task in Pro if you used a LYRX as the "Symbology Layer". But when I tried this it did not keep the break values constant (which is what I ultimately want).

In other words, if the layer I'm basing the symbology off of looks like this

I want the colors associated with 1-14 to remain constant.

Am I missing something? Is there a solution for this using ArcPy? Or is it just not possible unless ESRI adds the "Complete Symbology Definition" feature into ArcGIS Pro?


texts: How to reimport a QGIS 3.10 Vector-GeoPDF, generated from layout manager, so that all items not belonging to layers are hidden?

texts:
I'm very impressed with the new GeoPDF export capabilities of QGIS 3.10 (Windows 10). Now we can create Vector-GeoPDFs with layers included!

With the help of GDALs "Geospatial PDF" file format documentation (https://gdal.org/drivers/raster/pdf.html) I managed to reimport a GeoPDF, which I exported from the Layout Manager, turning on and off certain layers. Therefor I use DOS Environment variables to turn on (GDAL_PDF_LAYERS) and turn off layers (GDAL_PDF_LAYERS_OFF).

The GDAL command gdalinfo -mdd "all" helps to find all available layer names.

Unfortunately, I cannot hide all elements outside my layers (i.e. scale bars, legends, titles, etc.).

So my question:
Is there any hidden layername, or is there any GDAL trick to filter these items?

P.S.: I don't want to crop the image using the PDF NEATLINE entry, which seems to be incorrect anyway (the vertex coordinates of the polygon are in wrong order)


texts: Look up OS grid letters from easting/northing in R

texts:
Suppose I have an easting/northing point:

data<-c("423851", "300384")

I would like to find out the grid letters for this point using the British National Grid system (Ordnance Survey). If you visit https://gridreferencefinder.com/ and put those coordinates into the 'x (easting)' and 'y (northing)' boxes, it will tell you the two-letter grid - in this case SK.

Is there a package/function within r that will allow me to pass multiple coordinates and return the grid letters? I've searched for some time now and found nothing. I suspect that I'm just using the wrong terms but would like to check for a pre-existing function before considering writing my own.


texts: Issues with geometry in QGIS but not ArcGIS?

texts:
I have found that shapefiles that function normally in ArcGIS will not work in QGIS due to geometry issues. For example, I have a species geographic range file that when I intersect with a buffer will not work citing geometry issues (GEOS geoprocessing error. One or more input features have invalid geometry). When I try to fix the files in QGIS with various methods, Vector -> Geometry -> Check Geometries, convert to raster, create a buffer, manually edit out the self intersection I am unable to complete any of the processes because of geometry issues, or in the case of Check Geometries tool it runs to completion but does not seem to fix the problem. I have also run the GRASS v.clean and tried to manually edit out the self-intersection on a shapefile in QGIS but I am unable to edit the file. If I run any of these tools on the layer in ArcGIS it works no problem. This is my current solution, fix it in arc and then import it to QGIS.

My pointed question is therefore (1) have anyone experienced sensitivities to geometry in QGIS but not arcGIS, and (2) how do you fix files that seem to be rendered entirely unworkable in QGIS because of geometry issues?


texts:
Somehow, my attempt at setting the QGIS project coordinate reference system to a Sinusoidal projection by defining the CRS' EPSG ID returns Invalid CRS!.
crs = 54008 # World Sinusoidal Projection

# Set project coordinate reference system (CRS)
spatRef = QgsCoordinateReferenceSystem(crs, QgsCoordinateReferenceSystem.EpsgCrsId)

if spatRef.isValid():
    print("CRS Description: {}".format(spatRef.description()))
    print("CRS PROJ text: {}".format(spatRef.toProj4()))
    QgsProject.instance().setCrs(spatRef)
else:
    print("Invalid CRS!")

I use the long-term QGIS version 3.10.3 on Windows.


texts: Performing basic raster to shapefile conversion?

texts:
I am a novice when it comes to GIS. I am currently trying to use ERDAS imagine to classify habitat types using landsat 8 images. I'm able to run a supervised classification with the product being a TIFF file. I then opened this tiff file in ArcMap. Now I'm having a hard time projected this imagery over a base map in ArcMap or even extracting coordinates from the imagery with the respective classified habitat.

I tried to convert the raster to a vector/shapefile in Erdas but was given an error saying that the input file had to be thematic.


texts:
It seems to be that I'm too stupid to use gdal_grid. I have a
file called dgm-10-epsg-32633.csv with equidistant samples 10x10m in this form.

  X,Y,Z
  386700,6068670.000,.000
  386710,6068670.000,.000
  386720,6068670.000,.000
  386730,6068670.000,.000
  386740,6068670.000,.000
  386750,6068670.000,.000
  386760,6068670.000,.000
  ...

Before I start the computation I've calculated the extent.

Extent3d default {
  nanZ 9.999000e+03
  minX 3.700000e+05
  maxX 4.200000e+05
  minY 6.000000e+06
  maxY 6.068670e+06
  minZ -4.148000e+01
  maxZ 1.596600e+02
  count 31307010
  centerX 3.948742e+05
  centerY 6.031587e+06
  centerZ 6.993280e+00
  rangeX 5.000000e+04
  rangeY 6.867000e+04
  rangeZ 2.011400e+02
}

to get the parameters for -txe and -tye and -outsize.

I wrote a vrt file and the content looks like this

<OGRVRTDataSource>
    <OGRVRTLayer name="dgm-10-epsg-32633">
        <SrcDataSource>dgm-10-epsg-32633.csv</SrcDataSource>
    <LayerSRS>epsg:32633</LayerSRS>
    <GeometryType>wkbPoint</GeometryType>
    <GeometryField encoding="PointFromColumns" x="X" y="Y" z="Z"/>
   </OGRVRTLayer>
</OGRVRTDataSource>

and ogrifno gives me a valid response to data description.

$ ogrinfo -so dgm-10-epsg-32633.vrt
INFO: Open of `dgm-10-epsg-32633.vrt'
      using driver `VRT' successful.
1: dgm-10-epsg-32633 (Point)

If I call the commandline:

$ gdal_grid \
  -a nearest:radius1=10:radius2=10:angle=0.0:nodata=9999\
  -txe 3.7e+05 4.2e+05 \
  -tye 6e+06 6.06e+06 \
  -outsize 5000 6867 \
  -of GTiff -ot Float32 \
  -zfield Z \
  dgm-10-epsg-32633.vrt \
  dgm-10-epsg-32633.tif

I get always the help message without any explanation:

Usage: gdal_grid [--help-general] [--formats]
    [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
          CInt16/CInt32/CFloat32/CFloat64}]
    [-of format] [-co "NAME=VALUE"]
    [-zfield field_name]
    [-a_srs srs_def] [-spat xmin ymin xmax ymax]
    [-clipsrc <xmin ymin xmax ymax>|WKT|datasource|spat_extent]
    [-clipsrcsql sql_statement] [-clipsrclayer layer]

Available algorithms and parameters with their's defaults:
    Inverse distance to a power (default)
        invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0
    Moving average
        average:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0
    Nearest neighbor
        nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0
    Various data metrics
        <metric name>:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0
        possible metrics are:
            minimum
            maximum
            range
            count
            average_distance
            average_distance_pts

What is wrong on this stuff?


texts: ArcMap AddIn Build Error when referencing external ActiveX Interface in VS2010

texts:
If the "Register for COM interop" is unchecked the projects builds without errors.

But I have to reference external ActiveX Interface in my project so I have to have that check-box checked. I am getting the build error:

Cannot register assembly "...\MyCustomArcMapAddIn.dll". Could not load file or assembly 'ESRI.ArcGIS.Desktop.AddIns, Version=10.1.0.0, Culture=neutral, PublicKeyToken=8fc3cc631e44ad86' or one of its dependencies. The system cannot find the file specified.

Is there workaround here?
Thanks


texts: Writing standalone python script that opens canvas with map tools in QGIS?

texts:
I'm trying to write a stand alone script that will open a canvas with map tools and zoom to a lat/long with raster images. I'm starting by trying to run the script here http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/canvas.html#using-map-tools-with-canvas, and combining it with the standalone script here: http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/intro.html#using-pyqgis-in-standalone-scripts.

from qgis.gui import *
from PyQt4.QtGui import QAction, QMainWindow
from PyQt4.QtCore import SIGNAL, Qt

QgsApplication.setPrefixPath(""C:/OSGeo4W/apps/qgis"", True)

# create a reference to the QgsApplication, setting the
# second argument to False disables the GUI
qgs = QgsApplication([], False)

# load providers
qgs.initQgis()

class MyWnd(QMainWindow):
  def __init__(self, layer):
    QMainWindow.__init__(self)

    self.canvas = QgsMapCanvas()
    self.canvas.setCanvasColor(Qt.white)

    self.canvas.setExtent(layer.extent())
    self.canvas.setLayerSet([QgsMapCanvasLayer(layer)])

    self.setCentralWidget(self.canvas)

    actionZoomIn = QAction(QString("Zoom in"), self)
    actionZoomOut = QAction(QString("Zoom out"), self)
    actionPan = QAction(QString("Pan"), self)

    actionZoomIn.setCheckable(True)
    actionZoomOut.setCheckable(True)
    actionPan.setCheckable(True)

    self.connect(actionZoomIn, SIGNAL("triggered()"), self.zoomIn)
    self.connect(actionZoomOut, SIGNAL("triggered()"), self.zoomOut)
    self.connect(actionPan, SIGNAL("triggered()"), self.pan)

    self.toolbar = self.addToolBar("Canvas actions")
    self.toolbar.addAction(actionZoomIn)
    self.toolbar.addAction(actionZoomOut)
    self.toolbar.addAction(actionPan)

    # create the map tools
    self.toolPan = QgsMapToolPan(self.canvas)
    self.toolPan.setAction(actionPan)
    self.toolZoomIn = QgsMapToolZoom(self.canvas, False) # false = in
    self.toolZoomIn.setAction(actionZoomIn)
    self.toolZoomOut = QgsMapToolZoom(self.canvas, True) # true = out
    self.toolZoomOut.setAction(actionZoomOut)

    self.pan()

  def zoomIn(self):
    self.canvas.setMapTool(self.toolZoomIn)

  def zoomOut(self):
    self.canvas.setMapTool(self.toolZoomOut)

  def pan(self):
    self.canvas.setMapTool(self.toolPan)

The console is returning:

C:\OSGeo4W\bin\python.exe C:/OSGeo4W/Pytest2/canvasTools.py
ImportError: No module named site

Process finished with exit code 1


texts: Trouble with projections/inverse projections

texts:
I'm afraid this might be a very basic question but after scouring Stackexchange for a few hours and trying different things every few months I'm still out of luck. I'm trying to set a CRS / project my layers and with everything I've tried it is never straightforward and tells me (over and over) that this is an inverse transformation. I don't know how to project my layers, basically, and I've watched videos and followed lots of different articles about how and somehow it always ends up as an inverse transformation.
It seems like it should be a easy fix (aka, my ignorance on projections and coordinate systems is at fault) but I don't have anyone to ask about this so it's not the easiest to problem solve. Always after these inverse transformations, it says ballpark transform occurred. I'm using the correct coordinate system, pretty sure - my shapefiles should definitely be in zone 10N.
Things I've tried:

Save layer and change CRS in the save: EPSG:4326-WGS 84 to EPSG:26710-NAD27/UTM zone 10N [or NAD83/UTM zone 10N]
Go to layer>properties, try to change from there
Go to layer>set crs of layers
On layer > set CRS
project>project properties

An unrelated but maybe related problem - everytime I load QGIS I get an error message:
2021-02-08T16:13:08     WARNING    Python error : Couldn't load plugin 'processing' See message log (Python Error) for more details.

This is also something I've ignored for quite some time. Could this be related/causing my problems or is it a totally different thing?
(1)


texts:
I have a bunch of tracks on screen in a vector layer. I used some hover and select styles in the layer styleMap and these styles are properly triggered by:

    var highlightCtrl = new OpenLayers.Control.SelectFeature(tbTrack.vector, {
        hover: true,
        highlightOnly: true,
        renderIntent: "hover",
    });

    var selectCtrl = new OpenLayers.Control.SelectFeature(tbTrack.vector,
        {clickout: true, renderIntent: "select"}
    );

    tbMap.map.addControl(highlightCtrl);
    tbMap.map.addControl(selectCtrl);

    highlightCtrl.activate();
    selectCtrl.activate()

All these tracks are listed in a table under the map. I want to trigger the corresponding hover and select styles for a track if I click or hover on them in the table. How to do this?
I allready have the coupling between the table and the lines, I just don't know the right way to trigger the hover or select events.


texts: Field sampling design: stratified random or transect?

texts:
I am trying to understand the impact of scientific activity in a remote polar environment.

I have five locations in different parts of the region, where historic huts were placed for a couple of decades and then removed (approx. 1960s-1980s). Researchers would live in the huts and then go out and do surveys, and take measurements in the surrounding region. Each hut also has a designated site where helicopters would land and where and an outhouse was placed.

There are communities of microorganisms living in the soils, and these communities are understood to be highly vulnerable to perturbations.

My hypothesis is that due to human activity, soils in and around the footprint of the huts have more homogeneous biophysical characteristics than soils further away from the huts.

It seems there are two ways to approach this:

1) Stratified random Generating random sample points within the footprint and the area assumed to have human activity. This assumption would be based off of historic photographs and reports of the activities. A buffer would then be applied to create an outer ring and an equal number of samples would be randomly taken there. Each sample point could comprise multiple points in order to capture the variability of that point. Some kind of test would be run to compare differences between the in/out classes.

2) Transect Walk transects from starting from the center of the hut footprint or the centroid of the area assumed to have human activity. These could be three transects at three equally-spaced azimuths. Samples would be taken at regular intervals moving away from the hut area (again doing some kind of cluster at each one). This would allow me to get a sense of the gradient of change, but might make the calculation of a statistically significant difference difficult.

Please note that I have limited time and resources, so as few samples as necessary is ideal.

It also may not be possible to revisit sites, as these are highly protected areas, so adaptive sampling is likely not an option.


texts: Converting polylines to polygons in ArcGIS using multiprocessing

texts:
I am attempting to convert 15 feature classes containing roughly 70,000 contours each into polygons record by record. Using my current code, each feature class takes roughly 16 hours on my current computer using a Python script (i.e. not executed from within ArcGIS, only using the ArcPy library). My computer has a 3.2 GHz W3670 processer, 12 GB RAM, and an SSD. So, I think it should be able to process the data faster.

I have had no exposure to multiprocessing until I started looking around yesterday trying to find a method to speed up this process. I found several articles (see ESRI blog post and user post; earlier posts on GIS stack exchange 1 and 2)

Using the above examples, I attempted to mash-up what I saw there with my existing code and came up with:

import multiprocessing
import arcpy
arcpy.env.overwriteOutput = True

def contour2polygon(oid):
    contour = arcpy.management.MakeFeatureLayer('[some_path]\\test1.gdb\\contours','layer{0}'.format(oid[0]),'"OBJECTID" = {}'.format(oid[0]))
    out_name = "[some_path]\\test2.gdb\\feat{}".format(str(oid[0]))
    polygons = arcpy.FeatureToPolygon_management(contour, memory_name)
    return polygons

def main():
    contours = '[some_path]\\test1.gdb\\contours'
    oids = [row[0] for row in arcpy.da.SearchCursor(contours, 'OBJECTID')]

    cpuNum = multiprocessing.cpu_count() - 1

    pool = multiprocessing.Pool(processes = cpuNum)
    polygons = pool.map(contour2polygon, oids)
    pool.close()
    pool.join()

if __name__ == '__main__':
    main()

However, when I execute the script, I just see a blank box and my Task Manager does not indicate any change in activity. Is there something my code is missing? I understand from the above, linked examples that I should have a worker function, but I get lost in the rest of the details especially with how I could convert my old scripts to using multiprocessing.


texts: Finding and using data on population density files?

texts:
I'm attempting to replicate this map for the state of Nevada in QGIS:

source: http://www.bostonfed.org/commdev/c&b/2014/spring/map-distance-to-closest-snap-outlet.htm

(A SNAP outlet is just a store that accepts food stamps.)

I already got my zip code and census tract layers, and geocoded the SNAP outlets in Nevada.

Looking at the map, it doesn't seem like they used the centroid of each census tract to determine distance. I think they might have used population density. This makes sense because the majority of the population in any given tract could be say, in a corner. So it doesn't make sense to use the centroid of a tract. Assuming this is correct, what I want to do is import a population density point layer. Then, by census tract, find the centroid of the population density, and determine distance to the nearest SNAP outlet from there. I was able to find a population density grid, but I don't know how to use it. Upon importing it, it blackened my map.

I came across some population data on census.gov, but it appears to be population count by area (polygon), not point data.


texts: Coding names with a uniqueid starting form 001 but not exceeding three digits

texts:
I'm trying to code the names from geo database the total records of names are 957 and i want each name to have the uniqueID and it should not exceed three digits. starting from 001.

The below script started from 101 which now is giving me a problem when i rich to one thousand because it gives me the four digits.

I have been using the below Script to code names but unfortunately i'm now stuck because some ID'S are now having four digits.

rec=0
def autoIncrement():
 global rec
 pStart =101
 pInterval = 1
 if (rec == 0):
  rec = pStart
 else:
  rec += pInterval
 return rec

Here is the steps to code,

Right-Click on Field Calculator
Set the Parser to Python (Default is usually VB)
Tick Show Codeblock.

Paste the following into the Pre-Logic Script Code

rec=0
def autoIncrement():
 global rec
 pStart =101
 pInterval = 1
 if (rec == 0):
  rec = pStart
 else:
  rec += pInterval
 return rec

Paste the following code in the smaller box below the Pre-Logic Script Code:

autoIncrement()

Click OK.


texts: Finding Final Point using Azimuth and Distance Plugin

texts:
I am new to QGIS, only having used Arc before. I am working for a non-for -profit that is trying to track birds that are viewed from points hiking trails. The GIS intern before me used data on the bearing and distance of birds, and used the "Azimuth and Distance" plug in to create lines pointing to where the birds were. My boss wants these lines turned into single points where the birds were actually located. However, I do not know how to find the terminal point of these lines (the bird's real location) precisely.

So far, I have turned on the nodes of these lines, and found the terminal point, but this forces me to delete all the other nodes that also appear. Is there a plug in I can use to get precise end points for these birds?


texts: LAS dataset to raster with ground points only

texts:
I'm trying to make a raster from an LAS dataset. I want to only use points classified as ground, but I can't figure out how to do this. What I've tried is:

Create LAS dataset from 2 .las files.
Filter dataset on ground using LAS Dataset to raster tool.
Pick the dataset and use linear interpolation on average cell assignment type.

The result is a raster with all points including non-ground points.

When I'm in the LAS Dataset to raster tool I'm choosing an external file that doesn't get affected from what I've done in the LAS Dataset toolbar. How do I make the filter affect the dataset before turning it into a raster?


texts: ArcGIS server 10.3 fails to add machines to the cluster “Unable to access the config store”

texts:
Trying to Apply cluster scenario for ArcGIS Server on a test environment, 2 servers Node-02 and Node-03 these 2 servers are connected to a domain controller and have administrative users on the domain.

Now I am trying to add the machine Node-02 to the site that was created on the machine node-03, but I got the following error.

“Failed to register the server machine 'NODE-02.CLUSTER-MOLG.PS'. Server machine 'http://NODE-02.CLUSTER-MOLG.PS:6080/arcgis/admin' returned an error. 'Unable to access the config store on 'C:\arcgisserver\config-store'.'”

I believe that there is something wrong with ArcGIS server 10.3, since I can navigate to the folder 'C:\arcgisserver\config-store which is located on the machine Node-03 form the machine node-02, and the whole arcgisserver folder is shared with everyone with read and write property, so its not a network error.

What might be the issue here?


texts: Exporting a timeslice from a Netcdf-Rasterlayer as tif-rasterlayer

texts:
In ArcGIS Desktop 10.4 (advanced licence) I downloaded a NetCDF-File from here: http://eca.knmi.nl/download/ensembles/downloadchunks.php (the 0.25 grid with the years 1995-2016). It consists of temperature-data over time (in days).

I read it into ArcGIS as a NetCDF-Rasterlayer, and now I'm trying to export data of a single time-slice as a tif-rasterlayer. I did not get the hang on Mosaic Datasets, so my plan is to make a spatial analysis on the exported tif-rasterlayer, relating it to landuse (also in a tif-rasterlayer), and later aggregating it to NUTS3-Level). If I keep track what happens to which rasters I can later export other time-slices as tables and relate the spatial analysis to them.

The export-function (Data->Export) works in principle, but the raster looses it's attribute table, and without that I can't keep track which cell relates to which.

Does anyone have a suggestion on how to get that working?


texts: gdal_translate consuming too much memory when writing big JPEG2000 file

texts:
I am converting a 30.000 x 30.000 vrt image (made of JPEG2000 subimages) into a single (big) JPEG2000 image, using the KDU plugin.

gdal_translate -of JP2KAK -co QUALITY=100 my30000_3000.vrt myOutput.jp2

Now, my input is a 30.000 x 30.000 image, 4 bands, 8 bits per band. I would thus expect an average of 7GBi of RAM to be enough (let's add 1 GBi for "other processes"), but I would like to know if the conversion is possible with a limited amount of RAM, for instance, 4GBi.

I am trying to generate a 50k x 50k image with this method, but gdal_translate crashes after some time spent swapping bytes, keeping in mind I only have a limited amount of memory.

Is there a way to do so?

[EDIT] Versions : I am using gdal 2.2.2, with Kakadu 7.7[.1]


texts: Orient OpenLayers 6 labels along longest side of a polygon

texts:
I'm using OpenLayers with the EPSG 3857 Coordinate System and I want my polygon labels to be orientated along the longest side of the polygon.
The labels in OpenLayers are aligned horizontally by default but I can rotate them by passing a radian value. It goes clockwise on positive and counterclockwise on negative radian.
So I need to calculate the correct radians.
First I iterate the points of a polygon and calculate the longest side. Then I determine my vector from one point to another by subtracting the x and y values and the unit vector in the direction of the x-axis.
Now I calculate the scalar product and if it's zero then I know I have to rotate my labels by 90 degree or PI/2 for radian. This works fine I just add PI to flip the text so I can read it easier.
With if/else at the end I'm messing around to get the right angle but I'm not sure at this point because it works for the most polygons but not all of them.
What I'm missing here?
My math looks like this for each polygon:
var arrPoints = feature.getGeometry().getFlatCoordinates();

if (arrPoints.length >= 6) {

    var x = 0;
    var y = 1;

    var longestSide = 0;
    var longestSidePointFrom;
    var longestSidePointTo;

    for (var i = 4; i <= arrPoints.length; i += 2) {
        var pointFrom = [arrPoints[i - 4], arrPoints[i - 3]];
        var pointTo = [arrPoints[i - 2], arrPoints[i - 1]];

        const xDiff = arrPoints[i - 4] - arrPoints[i - 2];
        const yDiff = arrPoints[i - 3] - arrPoints[i - 1];

        var distance = Math.sqrt(xDiff * xDiff + yDiff * yDiff);

        if (distance > longestSide) {
            longestSide = distance;

            longestSidePointFrom = pointFrom;
            longestSidePointTo = pointTo;
        }
    }

    var vector = [longestSidePointFrom[x] - longestSidePointTo[x], longestSidePointFrom[y] - longestSidePointTo[y]];

    var scalarProduct = vector[x] * 1 + vector[y] * 0;

    if (scalarProduct == 0) {
        return Math.PI + (Math.PI / 2);
    } else {
        var sumA = Math.pow(vector[x], 2) + Math.pow(vector[y], 2);
        var sumB = Math.pow(1, 2) + Math.pow(0, 2);

        var cosine = scalarProduct / (Math.sqrt(sumA) * Math.sqrt(sumB));

        var radian = Math.acos(cosine);

        if (radian > (Math.PI / 2)) {
            radian = radian * -1 + (Math.PI)
        } else {
            radian = radian + 2 * Math.PI;
        }

        return radian;
    }
}

EDIT:
Here are 4 Polygons. The green line indicates the longest side. As you can see the angle for number 1 is not correct:

Here are the coordinates for the shown polygons:
1)
0: 44.7879
1: 313.8661
2: 42.1537
3: 317.1012
4: 45.4705
5: 319.8012
6: 48.1048
7: 316.5661
8: 44.7879
9: 313.8661
2)
0: 42.458
1: 311.9045
2: 39.7923
3: 315.1779
4: 42.1537
5: 317.1012
6: 44.7879
7: 313.8661
8: 42.458
9: 311.9045
3)
0: 37.2561
1: 307.6731
2: 34.5923
3: 310.9445
4: 36.9185
5: 312.8382
6: 39.5832
7: 309.5648
8: 37.2561
9: 307.6731
4)
0: 29.4927
1: 301.3244
2: 26.8156
3: 304.6124
4: 30.0662
5: 307.259
6: 32.7309
7: 303.9856
8: 29.4927
9: 301.3244
If I remove my if statements then it looks like this:


texts: Open KML Photos In Default Desktop Application

texts:
I have created a KML file containing several thousand placemarks.  The description contains a simple HTML  tag referencing a location on local disk, pointing to a JPG photograph.  The photos display in the popup balloons properly.

I would like for users of Google Earth 7 to be able to click on the displayed photo and have their default JPG file type handler application open and display the photo file.

I have tried enclosing the  tag with an  tag, using the notation href="file://photos\technician01\photo1.jpg", but I can't get the behaviour I am looking for (no external viewer is ever spawned).

Is there a way to accomplish this?

Update
Here's a complete block of HTML I have attempted to achieve the described functionality:

`


texts: gdal_retile is currently using only one CPU of the available eight and it is taking a while (48 hours already) to retile 5TB of data.


texts: I received a max limit on Geocoding address using MMQGIS when choosing OSM/Nominatim

texts:
I was trying to run a geocode process (~20,000 addresses) using MMQGIS and the OSM/Nominatim web service.

On the first run I input my CSV file and designated all parameters. However, I hit a max limit error for Google API.

Even though, I didn't choose Google API and my understanding of OSM/Nominatim  there are no limit to geocoding addresses.

When I tried to rerun the geocode service thinking it was a false error I got this.

I've tried searching for responses that may provide some insight, however I have been unsuccessful in finding anything useful. As far as I can tell, there should be no limit on the OSM and for Google is 2,500 which is fine, I'm not trying to use Google Map API web service.

Why am I receiving this error now?


texts: Documenting workflow and processing in QGIS-project

texts:
I want to include a documentation about the steps I undertook to create a QGIS-project to be able to understand how my data was prepared and processed, which steps I performed in QGIS. I'm not sure about the right place where this should be done. Looking to the Project Properties Dialog, section Metadata Properties, there is a history tab. However, the documentation does not really explain it:
https://docs.qgis.org/3.10/en/docs/user_manual/introduction/qgis_configuration.html#project-properties
Is this tab ment to get a description of the workflow or is there another place for it in a QGIS project?


texts: Leaflet mouse wheel zoom only after click on map

texts:
I'm working with the Leaflet JavaScript Library and attached a (working) map to my HTML Document. It is in the middle of the page, and when I'm scrolling down with my mouse wheel and arrive at the map, it automatically zooms into the map.

I want to scroll through the page without stopping at the map. Is there a way to activate the wheel zoom only after the first click on the map?


texts: QField App crashes when connect to internal server of PostgreSQL and PostGIS

texts:
I'm completely new to GIS field and here what I'm trying to do:

I want to make the QField App connects to internal server which has a PostgreSQL and PostGIS running on it.

I created a schema in PostgreSQL with one geometry type.
I enabled the PostGIS to the database
I created a QGIS project, and added the OpenStreetMap as a layer.
I also used Add layer and chose PostGIS and select my schema (from PostgreSQL)
Also, tried to export the project with QfieldSync

Also, I changed the IPv4 in pg_hba.conf and set firewall rules for PostgreSQL port(5432).

When I use QGis to select points, everything works fine and all data stored to my database. However, when I copy the project to my android phone and use QField, it crashes once I select the project without showing any error.

Is there a special way to export the project to fit with QField?

If I use shapeFiles as layer, the QField works just fine and store data locally but it is not what I need.


texts: Google map API with JSON response some point need to decoded

texts:
When I use Google map API with python, I  would like to know the path between the two points.
it response a JSON, but point need to be decoded.

example:

this is a set of point data : cq{oCqahzQDGr@iA

but how can I to decode this, making it become readable point coordinates?


texts: How to create different maps with same set of layers?

texts:
In ArcMap 10.7, I have built a map and created more than 40 layers e.g. for grassland, forest, water, borders etc.

Now I'd like to save that map as a template without any polygons or lines but with all the layers and symbols, because I need to create hundreds of different maps that all use the same layers.

So what I want is: Open an EMPTY template with the layer structure ready to go, import a raster img and start digitize the environment, finish the map, save and move on to the next one.

PROBLEM: The layers always have the same data source, even when I re-import a layer package every time I start a new map. So when I finished my first map, open up my template to start working on the next one, it's not empty but all the features of the first map are in there. And if I delete all polygons/lines from the attribute table in map 2, I also delete them from my first map because the reference is the same.

What am I missing??
I just need an empty map with all my (empty) layers loaded as a starting point.

One possible solution: I would always have to change the data source for every layer to a new shp-file, so the previous project remains intact. As soon as I forget to re-set the data source I would destroy other maps and hours of work.
Would also be a pain in the ### to change the data source every time.

PLEASE help. Even polyplot for DOS from the 90s can do it!


texts: Forcing v.buffer distance units to be different from map proj units using GRASS GIS

texts:
I'm not very familiar with GRASS GIS. I am working on MODIS data, and set the region to match the projection of the data, which is given in degree units. I need to create a buffer with a 250m distance, I tried to use as distance field 0:00:08.111518 (that is equal to the map resolution, that is 250m in my case), but it returns an empty map. I don't understand where is the mistake that I am doing,
What am I doing wrong?
I'd change the map units, if it can be done, but I don't know if it could mess the data up, therefore I was wondering if there is any method to force the unit of v.buffer to differ from the one of the map.


texts: Seeking New York Neighborhoods Boundaries?

texts:
I need a source of New York neighborhoods boundaries coordinates, I need to draw polygons around them just like this page is doing to Wingate, Brooklyn.
https://www.exrny.com/listings/121719/2-bed-in-east-flatbush-ny

I spent days looking on Google, the best I could find was api.maponics.com/sample, but they just give away ZIP boundaries on the sample and zip codes sometimes include more than one neighborhood.


texts: How to insert QGIS project properties "Project Title" as field or token using an expression?

texts:
Using QGIS, the "General" project properties dialogue includes a "Project Title" area.

Is it possible to use the information inserted in this dialogue in an expression (specifically including within a label when using Composer)?


texts: Enhanced Vegetation Index( EVI) Value Interpretation

texts:
Im working on calculation of Calculation of Enhanced Vegetation Index(EVI) From Landsat 5 TM and Landsat 8 OLI Using ArcMap 10.4.1. As far i know to Calculate EVI, first it needs to Convert the DN Values to Physical TOA reflectance. For Landsat 5 TM, I did that steps (First I converted DN to Radiance and then to Reflectance and  got the refelectance Values of the required bands ranging between  0 to 1 (Band 1 =0 to 0.367, Band 3= 0 to 0.648, Band 4= 0 to 0.808).

Then I applied the Formula for Deriving EVI= 2.5 * (("Band 4" - "Band 3") / "Band 4" + 6 * "Band 3" - 7.5 * "Band 1" + 1) in Raster Calculator and got result with EVI Ranging -1822.42 to 12.71 (!) . If i clip my Area of Interest Still get the Value ranges between -0.658 to +2.95.

I used to know that EVI ranges between -1 to +1 , However USGS says The valid Range is between -10,000 to +10,000) see https://www.usgs.gov/land-resources/nli/landsat/landsat-enhanced-vegetation-index?qt-science_support_page_related_con=0#qt-science_support_page_related_con
And i know EVI is  not a Normalized Index like NDVI.

Now I really Want to know if my precedure is okay, and then how can I interpret The values? (eg. Vegetation, Bare Soil etc).
And If not then how can I Correctly Calculate EVI?


texts: How do you add a search box to the map when creating a CartoDB visualization at runtime?

texts: I know that it is possible to add the CartoDB search box when you link to a visualization such as this. However, I have a map which is created at runtime, using the same method as here. Does anybody know, when a visualization is created at runtime, if there is any way to add the CartoDB search box?


texts: Possible to set an attribute field to update automatically based on another field?

texts:
If I have 2 attribute columns and I wish to change certain values, is it possible for the second attribute column to automatically update its values based on an expression?

Example:

Age |  Status
10  |  Child
14  |  Child

If I change the Age from 14 to 18, I would like the Status to automatically change from Child to Adult, so it would look like the following:

Age |  Status
10  |  Child
18  |  Adult

I must only be able to do this within QGIS and without the use of database functionalities such as PostGIS/SpatiaLite etc. I've looked at Value relations but that only seems to provide drop-down box for the user to select the relating columns. I've also looked at Relation Reference but it won't allow me to select a "Relation" (not even sure what Relation Reference actually does).

Please advise?

I am using QGIS 2.2.


texts: Selecting multiple attributes with variables using arcpy

texts:
I'm working on a program for my work and I'm trying to select multiple power line segments from the same shapefile. I'm having trouble with my SQL syntax to pass both variables on to the selection process in the program. I've tried both the 'or' and 'in' queries. The line of code worked with just one variable, but it still only passes one variable after I added the second one.  Any suggestions?

import arcpy

line_val1 = raw_input('Enter Line Number: ')
line_val2 = raw_input('Enter Line Number: ')

expression = str(""" "LINE_NUM_1" = '""" + line_val1 + "'" or """ "LINE_NUM_1" = '""" + line_val2 + "'")

arcpy.SelectLayerByAttribute_management(line_out, "NEW_SELECTION", expression)


texts: Generate Two Triangles on KML file with PHP/MySQL

texts:
I am currently generating a KML file dynamically from PHP and Mysql. One latitude / longitude point is being passed to the script. What I want to do now is create two triangle polygons from this point based upon a wind variable in the database. So if the wind is blowing 'S' then I want the triangle 'point' to start at the lat/long point variable, and build a polygon with the point of the triangle headed south towards the lat / lng. Can anyone share with me how this can be done? Here is my script:

// Selects all the rows in the markers table.
$query = 'SELECT lat, lng, waves, wave_secs, current_dir, pri_id AS id FROM kml';
$result = mysql_query($query);

if (!$result)
{
die('Invalid query: ' . mysql_error());
}

// Creates the Document.
$dom = new DOMDocument('1.0', 'UTF-8');

// Creates the root KML element and appends it to the root document.
$node = $dom->createElementNS('http://earth.google.com/kml/2.1', 'kml');
$parNode = $dom->appendChild($node);

// Creates a KML Document element and append it to the KML element.
$dnode = $dom->createElement('Document');
$docNode = $parNode->appendChild($dnode);

// Iterates through the MySQL results, creating one Placemark for each row.
while ($row = @mysql_fetch_assoc($result))
{
// Creates a Placemark and append it to the Document.

$node = $dom->createElement('Placemark');
$placeNode = $docNode->appendChild($node);

// Creates an id attribute and assign it the value of id column.
$placeNode->setAttribute('id', 'placemark' . $row['id']);

// Create name, and description elements and assigns them the values of the name and address columns from the results.
$nameNode = $dom->createElement('name', 'Wave Height '.htmlentities(ceil($row['waves'] * 3.280)).' Feet');
$placeNode->appendChild($nameNode);
$descNode = $dom->createElement('description', 'Wave Duration ' .htmlentities($row['wave_secs']). ' Seconds');
$placeNode->appendChild($descNode);
$styleUrl = $dom->createElement('styleUrl', '#' . $row['current_dir'] . '-Style');
$placeNode->appendChild($styleUrl);

// Creates a Point element.
$pointNode = $dom->createElement('Point');
$placeNode->appendChild($pointNode);

// Creates a coordinates element and gives it the value of the lng and lat columns from the results.
$coorStr = $row['lng'] . ','  . $row['lat'];
$coorNode = $dom->createElement('coordinates', $coorStr);
$pointNode->appendChild($coorNode);
}

$kmlOutput = $dom->saveXML();
header('Content-type: application/vnd.google-earth.kml+xml');
echo $kmlOutput;


texts: Stroked and filled vector lines in Openlayers

texts:
I am drawing vector lines, using GeoJSON, in my web application that uses Openlayers. I want the line to be displayed in yellow color with black border. I am using the following for the line styling:

my_layer = new OpenLayers.Layer.Vector("my_layer", {
            styleMap: new OpenLayers.StyleMap(new OpenLayers.Style({
                    fillColor: "#ffff00",
                    strokeWidth: 6,
                    strokeColor: "#000000"
            }))
        });

But 'fillcolor' does not affect line symbology. How can I achieve that effect?


texts: QGIS virtual layers and slow performance

texts:
I am using QGIS to create a layer comprised of counties in 7 northeast US states.
I start with the US census 5M cartographic state borders file, and then set a filter on that layer:
"NAME" in (
  'New York',
  'Massachusetts',
  'Rhode Island',
  'Connecticut',
  'New Hampshire',
  'Vermont',
  'Maine'
)

That gives me a layer with seven features. I wanted to take the US census 5M cartographic county borders file and restrict that to the same geographic area. I thought I could do it by setting up a virtual layer and setting up a spatial join, like this:
select a.*
  from cb_2018_us_county_5m as a
  join cb_2018_us_state_5m as b
  on (st_contains(b.geometry, a.geometry))

That works, but it takes 5 or 6 seconds to render the county outlines; I can watch as they show up one by one. Given the relatively small number of features involved in this operation, I was surprised by the performance. Am I doing something suboptimal? The same query in, say, PostGIS is a sub-second operation.


texts: Split polygon into four with respect to single point

texts:
Some days ago one of my colleague came to ask me a query related to Geometries.

 The situation is as.... There is one Polygon (suppose a city polygon) and there is one point (suppose user's current location) inside that polygon.

Now, what I have to achieve is to split polygon into four quadrant so as the lines which splits polygon into four should be perpendicular to each other and should pass from point. (please see image to help with understanding).

In other words, this processing will generate four section of polygon as NE,SE,SW and NW with respect to point(user's location).

So, how to achieve it? One immediate thought in my mind was to
1]Compute the bounds of polygon
2]Get distance by max-min/2
3]Use this as a radius of buffer which will drawn from point.
4]Then get the intersection of this buffer and original polygon
5]From these intersection to point draw a lines to split the polygon.

Unfortunately, it's need not be always right (simply bcoz point could be anywhere inside polygon).

It would be great if anyone find a way to solve this problem.

I'm trying to solve this problem using opensource Geospatial libraries and software's.


texts:
I am running GRASS GIS 7.0.4 on Ubuntu 16.04.

After completing a tutorial on extruding vector data, I tried running the nviz command to see the result in a 3D-view.

This is in norwegian, but it tells me that the folder/file nviz does not exist. What am I missing here?


texts:
I am adding a number of spatial csv files as layers and want to assign a qml layer style.

import processing

# User Dialog to pick Import_QGIS folder
Folder_Path = QFileDialog.getExistingDirectory(caption="Pick the folder 'Import_QGIS'.")

# Folder_Path for uri
Folder_Path_uri = "file:///"+ Folder_Path + "/"

# Make a tuple with the filenames
import os
i=0
my_tuple=[]
for file in os.listdir(Folder_Path):
    First3Letters = file[0:2]
    try:
        int(First3Letters) < 300
        print(file)
        my_tuple.append(file)
        i+=1
    except:
        continue
iNum = int(i)
print(iNum)

for x in range(iNum):
    try:
        uri = Folder_Path_uri + my_tuple[x] + "?type=csv&delimiter=,&detectTypes=yes&xField=X_Wert_Pos&yField=Y_Wert_Pos&crs=epsg:25832&spatialIndex=no&subsetIndex=no&watchFile=no"
        print(uri)
        mem_layer = QgsVectorLayer(uri, my_tuple[x], 'delimitedtext')
        QgsProject.instance().addMapLayer(mem_layer)
        mem_layer.loadNamedStyle("C:\Layer_styles\style3.qml")
        x+=1
    except:
        break

All layers in the specific folder are added as new layers, but only with the default style. What am I missing?


texts:
First of all I'm new to Tilestream, MBTiles and GIS.
As a part of of learning purpose I set up Tilestream with the help of tilestream tutorial and copied two MBTiles file to the tiles folder I am able to navigate to the tiles with the help of default Tilestream. http://localhost:8888
But when loading the tile url it says
Error: Tile does not exist
    at Statement. (/home/../tilestream/node_modules/mbtiles/lib/mbtiles.js:168:33)
can any one help me


texts: I am new to this. My Giswater plugin icon does not function upon clicking. But successfully installed. I have installed postgresql, pgadmin 3, enabled spatial extensions postgis and pgrouting,


texts: Error while installing GeoPandas libraries using prompt

texts:
I am getting this error when trying to install  GeoPandas libraries.....
Anyone with a solution?
(base) C:\Windows\system32>conda install -c conda-forge geopandas
Collecting package metadata (current_repodata.json): done
Solving environment: failed with initial frozen solve. Retrying with flexible solve.
Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source.
Collecting package metadata (repodata.json): done
Solving environment: failed with initial frozen solve. Retrying with flexible solve.
Solving environment: |


texts:
I would like to use grass7:i.image.mosaic algorithm in R using RQGIS packages But I don't understand how I can define it

I have 14 rasters named r1:r14. I do

params <- get_args_man(alg = "grass7:i.image.mosaic")
params$input <- paste0("r",1:length(my.raster.list))
params$output <- file.path(tempdir(), "ger_coords.tif")
params$GRASS_REGION_PARAMETER <- c("277743.101135","712479.268554","1590501.95667","1913668.6652")

out <- run_qgis(alg = "grass7:i.image.mosaic",
                params = params,
                load_output = TRUE)

But the GRASS_REGION_PARAMETER seems to be wrong

Error in run_qgis(alg = "grass7:i.image.mosaic", params = params, load_output = TRUE) :
  Error: Wrong parameter value: 277743.101135712479.2685541590501.956671913668.6652

I have taken a look at this answer but it is still dark. I need a hand ! I'm using R version 3.4.0, QGIS 2.18 and Grass 7.2


texts:
I'm trying to convert data from an AUtoCad Civil 3d drawing to a text format to use in ArcGIS. Each dwg has coordinate and text information that I would like to use for point location and attributes. The text looks to be a table but is in fact just mtext and line objects.
I have found LISP code that will export all the mtext to a text file but I am wondering if there is a way to place all the text found at the same xy location in the dwg to be placed on one line in the text file?

Thanks!


texts: In OpenLayers how do I prevent linestrings and polygons disappearing as I zoom out?

texts:
I can see how to do this using a feature style function with 2 styles (similar to the Custom Polygon Styles example). One would be the normal linestring/polygon one and used at high enough zooms.

The second would show just the first coordinate as a point and used at low zooms. I already use style functions, but none have depended on resolution yet.

The cut-off point between the two styles would depend on how big the feature would appear at that zoom (i.e. on resolution, and its extent?).

This seems a standard thing that anyone would want, so is there a simpler way of doing it? If not, can you give me any guidance on how to set the cut-off point?

I have seen the issue addressed elsewhere, but the only answer seemed to be about GeoServer, which I am not using.


texts: Converting text strings to numbers using QGIS?

texts:
I am working in QGIS, which is quite new to me.

I have a text string containing dates.

I want to have new column where only the year from the text-date string appears. How is that possible. See picture below:


texts:
I have glacier thinning data from https://github.com/dshean/hma_mb_paper/tree/master/data/mb in the form of CSV, which contains essentially everything, dh, dt all the parameters required.
I am just a newbie who'd want to create a map shown in the picture based on the data from CSV and RGI(Randolph Glacier Inventory that'll provide the glacier outlines in the form of polygons).
How do I do that? I am only aware about adding points (in this case, glaciers from x, y data) from the csv.


texts: Warning inconsistent extent when add shapefile


texts: Get displaylevels/zoom level javascript API

texts:
Is there any way using JavaScript API to get the current "displayLevels" of the basemap or "basemap zoom level" or "basemap scale" of the map and if it was more/less than certain number I change the Layer to another layer?

function(
        Map,
        FeatureLayer
        ){

map = new Map("map", {
           basemap: "topo",
           center: [-85.416, 49.781],
           zoom : 5,
           logo: false});
if (map zoom level == ***) {

 layer1 = new FeatureLayer(
                "http://testServer/arcgis/rest/services/test/MapServer/0", {
                        mode: FeatureLayer.MODE_ONDEMAND
                        });

        map.addLayer(layer1);
}

else {
layer2 = new FeatureLayer(
                "http://testServer/arcgis/rest/services/test/MapServer/1", {
                        mode: FeatureLayer.MODE_ONDEMAND
                        });

        map.addLayer(layer2);

}


texts:
I am building a small multispectral camera based on the Raspberry Pi Compute Module and two Raspberry Pi cameras. One of the cameras has a NIR filter, the other is a normal RGB camera.
When using the Composite Bands tool in ArcMap, and mix the bands as to get a falso color image, I got a blue image, rather than red. The combination is NIR, RED, GREEN.

That's the end image.

Any ideas why do I get this blue instead of red color?


texts: I'm using PostGIS 2.0.0 and would like to store a GML Surface in a way that allows calculating intersections between that surface and various points and lines. In my case, the Surface will have one PolygonPatch with no interior rings and an exterior ring which may be a LineStringSegment, a Geodesic, an Arc or a Circle. The distances can be in the order of tens to hundreds of miles. I see that PostGIS has CIRCULARSTRING and CURVEPOLYGON, which seem to be a good way of storing an Arc, except that they're not supported by the geography type, only geometry. What would be the best way to store this? (I'm quite new to PostGIS and GIS in general, so it's entirely possible I'm missing something obvious.)


texts:
My DEM raster coordinates are in WGS84: 45.4763884,-121.5736061 and vertical units are in feet I believe. I want to do a slope analysis.
Do I need to add scale factor?

I have been seeing different answers all over and am confused.


texts: Removing extra string from text column in Attribute Table to make numeric column from number part using ArcGIS Field Calculator?

texts: I received a shapefile about invasive species. I am using ArcGIS 10.2 for desktop. In Attribute Table, there is column called "Infested_Area" and its got number mix with text, for instance "9100 sq feet". Since I can not run statistic on string I was wondering if it's possible to create a new filed and somehow get rid of the text portion of it and only have number!


texts: ESRI Javascript API OAuth support for non-arcgis-online application

texts:
We are just in the opening salvos of developing a self-hosted ArcGIS Web application. All the OAuth examples I see seem to be assuming ArcGIS Online logins, or that the application will use those services.

But knowing how OAuth is supposed to work, we should be able to implement it in our IIS application and authenticate with Google, Facebook, or other OAuth providers.

How much of our own code will we have to roll for this to happen?

Am I missing something basic in the examples I've seen?

Can just snap together the existing pieces for our application without involving ESRI's servers, and without writing a bunch of C#?

Obviously implementing our access policies is on us. But I'm curious about the rest of it.


texts: Field Calcuator to add percent using Python

texts:
What might I go wrong here is that I am trying to add percent. I use the Field Calcualtor and it keeps getting an error.

I am opening the Table to add a field named Percent and then type this:

round([T_Acres]/207 * 100,0) + "%"


texts:
I downloaded a Digital Elevation Model of Mexico from here, and when I extracted the .rar, I got these files: Republica30_R60m.bil, Republica30_R60m.hdr, Republica30_R60m.prj and Republica30_R60m.stx.

I need to be able to read this for an Android App, but I haven't found any information about how to read these files.

Can anyone tell me how to read it or where I can find information about it?


texts: Using .tif on data layer for Google Maps?

texts:
i have three files that i should use for a GIS application: map.aux, map.tfw, map.tif

The .tif file is a "layer" that highlights some building over a map.

How can i bring this information (and "image") on a layer of Google Maps? Any idea?


texts: What's the most size-efficient spatial database

texts:
I'm configuring a tile server on Apache Mesos. From all the tutorials I've read, I gathered that PostgreSQL and PostGIS is widely used for this purpose.

I learned that MongoDb and MySQL Server also can be configured for spatial databases (but could not find any tutorial on the same).

My goal is to set up a tile server (which I've done with PostgreSQL) but wanted to know which DB would be the most efficient (w.r.t space and tile rendering speed)?


texts: Use Postgis to convert wkb_geometry data type to geom datatype

texts:
How would I write an ALTER statement or otherwise to convert an existing imported data column from "wkb_geometry geometry(LineString,4326)" to a geom data type?

I ask because I've moved to importing my data from geojson instead of shapefiles due to column length issues with shapefiles and now I've created the new problem for myself of the data type not being geom as default.


texts: When subtracting two raster with different cell sizes using raster math, what value is used for smaller cells?

texts: After reading the answer to a similar question (What will happen if we subtract two raster with different cell size using Raster Calculator?), I am still left with a question. When the raster with smaller cells is resampled, how do the values in the smaller cells get aggregated into a larger cell? Does it take the mean of the smaller cells being combined? Is there a way to choose?


texts: Tool for merging OSM and Google route data

texts:
I'm using google API in order to solve a routing problem; by giving it points of source and destination and demanding the route. Now I want to project the route on an OSM map.

Is there any tools to smooth the result and make sure that the route hasn't lost it's precision during the conversion?

I'm currently trying to check every segment of the route and if both ends lied on the same street, I show the segment.


texts: Maritime route prediction avoiding land masses

texts:
I am working on a project where the following functions have to be implemented:

Predicting the locations of the ships (in maritime environment) at a future time. (Can be done with Kalman filter, IMM filter and some other algorithms.) Ships can be in any part of the world.
Avoiding landmasses during prediction.
Finding the shortest paths along the shorelines.

I am totally done with the first part which is predicting without considering the shoreline information. I have problems with the functions 2 and 3.

Function 2 At times, the predicted location can fall into the landmass area which is totally unacceptable. I am using the coastal area shapefile http://openstreetmapdata.com/data/coastlines. This file has converted XY values of the world shoreline data. I have loaded this shapefile into postgreSQL and used postgis to read it from the database.

My idea is to go through all the polygons (shoreline defined based on polygons) and check whether the line connecting the present location and the predicted location crosses the polygon. If it crosses, that means we have to find where the ship intercepts the shoreline first. But if I follow this approach going through all the polygons, it is going to take forever. (It has around 62000 polygons and each of them has thousands of points). Any advice on this?

I thought about initially dividing the world map into hierarchical areas (Level 1: 10 polygons; Level 2: each polygon has 10 polygons inside). But I am not sure how to divide the world map into the levels of polygons I require. Is any functionality of postgis helpful for this? Or any other libraries for this purpose? I believe this kind of functionality should be available already. But I am not able to figure it out so far.

Function 3 Since now we know where the ship intercepts the shoreline first, we can predict its route along the shoreline using the shortest path algorithm given we know the destination information. But to do this, you need to divide the shoreline map into grids so the shortest path can be used. So how can you make grids based on this along the shorelines? I am not doing image processing here. What I have is this shapefile now. Or should I go with some image processing approach and make the grid shorelines? If so please provide some links.


texts: How to make a database query layer editable in QGIS

texts:
I have a layer in my QGIS 3.4 project that comes from a Geopackage database. I load it using a query to select only a subset of columns / fields and to filter partial records.

However, the layer cannot be edited on the canvas and the only way I can see to edit the data is to use dbManager and load the whole table (right click 'Add to canvas'). Unfortunately my table has many columns, not all of which are relevant to the 'canvas' / front end, so it makes it more difficult to work with the data.

Is it possible to edit the attributes of a layer that is loaded from a query?

I recognise that if I have multiple layers sourced from one table there is a risk of conflicts if I were to try editing all the layers concurrently... Maybe we can skip over that for now, if I promise to be careful!


texts: Visualizing lat/lon pairs on a map, locally

texts: I have a sensitive list of coordinates that I've already managed to geocode locally and now need to visualize on an interactive map. They're all contained within a particular neighborhood of the same city, so I'm hoping there's a way to display these on a map as pins, like one would do with the Google Maps API, but by downloading the map of the area and then adding the pins rather than sending the list of the pins and their data over the Internet. How might I get these points visualized on a map in a way that preserves the privacy of the coordinates?


texts:
As of yesterday I can no longer create new features in shapefiles that were previously working fine. There's no error messages, I can go through the process of adding a feature as normal but it just doesn't create it - in either the map or attribute table. Any solutions?

qgis v2.18.3


texts: Is it possible to dynamically construct the SQL query where clause in ArcGIS layer based on the URL parameters?

texts:
I have an ArcGIS layer published to my local server at:

https://baseURL/myMapServer/0

I want to dynamically construct the query where clause for this layer based on the URL. For example:

https://baseURL/myMapServer/0?category=cat1

would effectively translate into the following where clause:

category='cat1'

How can I accomplish this?

The arcgsis documentation touches on the topic:

https://pro.arcgis.com/en/pro-app/help/data/query-layers/update-parameter-values-in-a-query-layer.htm

https://www.esri.com/arcgis-blog/products/sharing-collaboration/sharing-collaboration/using-url-parameters-in-web-apps/

https://developers.arcgis.com/labs/arcgisonline/query-a-feature-layer/

And I've tried:

https://baseURL/myMapServer/0/?query=category=cat1

and then,in definition query:
category=::category

But it doesn't work, ArcMap throws an SQL syntax error.

Update: this question is specifically on how defining ArcGIS layers. This is not a question on how to retrieve GET query parameters in PHP, JavaScript or any common programming langauge.


texts: CSVLayer as a source for Search widget in ArcGIS Javascript

texts:
I am trying to add my CSVLayer as a source for the Search widget, however, I've been unsuccessful.  My Search widget searches addresses and I was able to add one of the sources from this example (https://developers.arcgis.com/javascript/jssamples/search_multiplesources.html) successfully, however my CSVLayer doesn't return any results in the Search widget.

Has anyone had any success or experience with adding a CSVLayer as a source for the Search widget in ArcGIS Javascript?  Below is a code snippet...

var csv = new CSVLayer("file.csv", {
      outFields: ["*"]
    });
var search = new Search({
      enableButtonMode: false,
      enableLabel: false,
      enableInfoWindow: true,
      showInfoWindowOnSelect: false,
      map: map
    }, "search");

    var sources = search.get("sources");

    sources.push({
      featureLayer: csv,
      searchFields: ["Name"],
      displayField: "Name",
      exactMatch: false,
      outFields: ["*"],
      name: "Foo",
      placeholder: "FooPlace",
      maxResults: 6,
      maxSuggestions: 6,

      infoTemplate: template,
      enableSuggestions: true,
      minCharacters: 0
    });

Could the issue be that my csv var is technically a CSVLayer and not a FeatureLayer?


texts: Collapse rows with different attributes in QGIS modeler

texts:
I am trying to build a QGIS modeler tool to do the following;

I have a parcel table with fields "OWNER" and "OWNERNO" where each owner has their own row in the table and the "OWNERNO" numbers them. (All other fields are the same in the rows.)  I would like to collapse this table so that:

Owner=Bill, OwnerNo=1
Owner=Mike, OwnerNo=2

Becomes:

Owner1=Bill, Owner2=Mike

I have searched this forum but could not find an answer.


texts: How to modify the GeoServer preview page - format=application/openlayers&layers

texts:
I want to modify the GeoServer preview page that brings up a simple OpenLayers page.

The GeoServer documentation is very rare: OpenLayers format=application/openlayers
Just saying: 'Generates an OpenLayers HTML application.'

Where and how can I modify this page? For example, I want to expand the map preview.


texts:
I am applying FLAASH to atmospherically correct and convert to reflectance a hyperspectral image acquired by our airborne system. The problem is that FLAASH (using ENVI) outputs an image where all spectra peaks around 763nm. all spectra have the same peaks regardless of the observed materiel. the Image bellow is a spectrum of a Spectralon that should have 36% reflectance.

My input image to FLAASH is Radiance acquired by the airborne system but actually has a weird problem too. Spectra regardless of the observed material drop in a negative peak at the same wavelength. Here is an image of the same Spectralon before applying FLAASH :

I am starting to think that this is a camera problem. Does anyone know where this problem is coming from? or suggest any ideas?

PS : I've applied FLAASH with different parameters but the peaks are always there.


texts:
When I use the osm2pgrouting I only get emtpy tables.

I got my .osm file from this Wien.os.gz and use the mapconfig.xml from Github.

I run the follwing command:

./osm2pgrouting -f path/Wien.osm -c path/conf.xml -d osm2 -h 127.0.0.1 -U postgres -W user --clean

I added the

 -h 127.0.0.1

according to this post, but it didn´t help.

This is the output I get:

Any ideas why this doesn´t work..


texts: Shapely intersection operation gives non-intersecting results


texts: Where are wrappers for .NET and ArcGIS 10x?

texts:
With prior versions (to ArcGIS 10x) of ArcGIS, the .NET wrappers used to be placed in a dotnet directory.

Where are they located with 10x?

With prior versions, I used the absence or presence of this directory to tell whether or not the .net framework was installed in the proper order with installation of ArcGIS.

How can I tell if the .net framework was installed in the correct order with respect to ArcGIS?


texts: Geographic Heatmaps with Concentric Circles by Distance?

texts:
I'm looking to generate a geographic heatmap (using 'ggmap') that overlays some dimension (to start, housing prices) over the lat/lon near a city center. Then I want to create circles of equi-distant spacing (i.e. 10 miles per circle) to get an idea how far out. I would also like my heatmap to go from blue to red for low to high of the dimension. I've been struggling with this for a day and this is as far as I got:

require(ggmap)
require(ggplot2)

ggmap(NewYork)
  + stat_density2d(data=positions, mapping=aes(x=lon, y=lat, fill=..level..), geom="polygon", alpha=0.2)
  + geom_point(shape=1, aes(x = housing.data.NY$Longitude, y = housing.data.NY$Latitude, size=sqrt(distance)), data = positions, alpha = .9, color="black")
  + scale_size(range=c(3,20))
  + labs(x = "Longitude", y = "Latitude", fill = "Housing \n Price Density")
  + ggtitle("Housing Prices by Distance\n New York, December 2014")

The code does the following:

Load the created GoogleMap file as a layer
Create price heat maps
Add concentric circles with radius ~ distance from city center (NEEDS WORK)
Scale the circles (or atleast try to)
Add labeling to make the plot more legible
Add plot title

The code produces the following output:


texts: Comparing two shapefiles(layers) and display the differences in features in MapContent object in GeoTools

texts:
I am working with GeoTools libraries to implement GIS features. Here, I have two different shapefiles and I want to generate a new separate layer containing the differences between the two shape files. This new layer should show the extra features which are there in one of the two input shp files.

It generates a new layer showing the differences between polygon2.shp and polygon1.shp.

Are there any specific methods in the GeoTools libraries which can capture the differences if we provide the two layers of shapefiles?


texts: I am designing a number of maps for a report and I was wondering is there a location where I can download different layout templates as the default ones in arcmap don't suit what I'm looking for?