texts: Finding the angle between two GIS linestring features using R

texts:
I have two simulated linestrings, one that represents a satellite ground-track and the other representing a linear surface feature (e.g., sea ice lead, road or river). I'm trying to determine how to calculate the inner angle of two line features. My eventual goal is to simulate many "lead" features and calculate their distance from and angle relative to the transect.
From searching stack exchange, I've been unable to identify a solution for this problem.
My current code approach is as follows:
library("sf")    # classes and functions for vector data
library(spData)        # load geographic data
library(lwgeom)
library(spatstat)
library(ggplot2)

#Simulate transect coordinates
x0 <- runif(1, 1,10)
y0 <- runif(1, 1,10)
x1 <- runif(1, 1,10)
y1 <- runif(1, 1,10)

#Simulate lead coordinates
l_x0 <- runif(1, 1,10)
l_y0 <- runif(1, 1,10)
l_x1 <- runif(1, 1,10)
l_y1 <- runif(1, 1,10)

# Set up a dataframe for the ICESAT-2 transect
transect.df <- data.frame(
  'x' = c(x0, x1),
  'y' = c(y0, y1)
)

# Set up a dataframe for the simulated lead info
lead.df <- data.frame(
  'x' = c(l_x0, l_x1),
  'y' = c(l_y0, l_y1)
)

transect.sf <- st_sf(id = 't01', st_sfc(st_linestring(as.matrix(transect.df)))) #create a formal sf feature
st_crs(transect.sf) <- st_crs(4326) # assign crs
transect.sf <- st_transform(transect.sf, crs = 3031) # transform

lead.sf <- st_sf(id = 'L01', st_sfc(st_linestring(as.matrix(lead.df)))) #create a formal sf feature
st_crs(lead.sf) <- st_crs(4326) # assign crs
lead.sf <- st_transform(lead.sf, crs = 3031) # transform

st_write(transect.sf,
         "transect1.shp",
         delete_layer = TRUE)

st_write(lead.sf,
         "lead1.shp",
         delete_layer = TRUE)

#Check distance
st_distance(transect.sf,lead.sf)

#Check angle
m1 <- ((y1 - y0)/(x1 - x0)) #transect slope
angle.m1 <- atan(m1)
m2 <- ((l_y1 - l_y0)/(l_x1 - l_x0)) #lead slope

angle <- atan(abs((m2-m1)/(1+m1*m2)))

ggplot() +
  geom_sf(data = transect.sf, fill = "grey", color = "red")+
  geom_sf(data = lead.sf, color = "blue")+
  ggtitle("Simulated Leads and ICESat-2 Transect") +
  coord_sf()


texts: After Merging Tiff Files the merged file is black and gray

texts:
I have 24 Tiff Rasters which I want to combine into one TIFF file but when I merge them the color of the new TIFF is black and gray.  The render type of the new TIFF is singleband gray but if I change it then it turns funky colors and does not look correct.

I merge by going to processing > toolbox > GDAL/OGR merge


texts: QGIS 3.4.x Properly digitizing polygons of uneven sides

texts:
I am facing some problems when digitizng polygons of uneven (irregular) sides using advanced digitizing and none/planimetric measurements for area and length calculations.

In the example image below the blue lines were created using advanced digitizing with locked distance (5 lines measuring 20.00 meters and one line 19.90 meters) and using perpendicular and parellel relative angles, the red ones were created also using advanced digitizing but they were the result of the intersection of two circles (20.00 meters, and 20.10 meters) using 'add circle by a center point and another point', both circles made at the end of the corresponding lines snapping on vertex, the green line was created snapping on vertex at the end of both lines. The length and area displayed is calculated by QGIS.

The workflow I just described is similar to what I do in AutoCAD, in it I can draw the same lines and polygons and get 400.00, 399.98, 400.00 square meters for P-1, P-2, P-3 respectively. I do not want to go back to AutoCAD because the print layout for QGIS is far superior to the one in AutoCAD (using atlas).

My problem is mostly with digitizing the lines using circles, the red lines should be of the proper length (20.10 meters and 20.00), so that way I can digitize P-1 again using parallel and perpendicular relative angles and get a proper value of 400.00 sqare meters.

The measurements for the polygons are taken from property deeds and titles (cadastre archive of local goverment), and the polygons are used to print parcel data that is given to the owner, so I cannot change the area of polygons that have even sides described in the property deeds, I could simply enter manual data for such polygons but I want to avoid that because I consider it malpractice. And this is why I do not want to go back to AutoCAD, as the printed document from AutoCAD requires manual entry of data (owner name, parcel code, zip code, etc.) and is prone to typing errors, the data is already available in a database that I can use for the polygons.

I also want to avoid having to import dxf data.

So my real question would be, what is the proper or better way of digitizing polygons of uneven sides using known segment lenghts?


texts:
I want to install the netCDF plugin for Geoserver, in order to publish radar data in netCDF format. I'm following the instructions listed at Geosolutions IT Training.
I have to add

-Duser.timezone=GMT

and

 -Dorg.geotools.shapefile.datetime=true

to the setenv.bat.

As stated in the training:

  Support for timestamp is not part of the DBF standard (used in
  shapefile for attributes). The DBF standard only supports Date, and
  only few applications understand it. As long as shapefiles are only
  used for GeoServer input that is not a problem, but the above setting
  will cause issues if you have WFS enabled and users also download
  shapefiles as GetFeature output: if the feature type extracted has
  timestamps, then the generated shapefile will have as well, making it
  difficult to use the generated shapefile in desktop applications. As a
  rule of thumb, if you also need WFS support it is advisable to use an
  external store (PostGIS, Oracle) instead of shapefile. Of course, if
  all that’s needed is a date, using shapefile as an index without the
  above property is fine as well.

I'm not sure I understand this correctly..so my question is:

Adding datetime=true in the setenv.bat will affect my capability to download WFS features in a .shp format, even if they don't have a Date field?


texts: Count classified attributes using second series of attributes in QGIS

texts:
I have:

classified GPS-data them according to "Starting city - Destination city".
a geodataset for German highways, that contains a "Route-ID" for each section.

Now I want to count how many times I drove over which "Route-ID" at which relation in QGIS.

For this, I joined the "Route-ID" per location and kept only one GPS-point per trip per "Route-ID". So I want to count the GPS-points per Relation per "Route-ID". Any suggestions?

The Dataset looks like:

As example:
I drove 5 times from Cologne to Frankfurt and used 4 times the combination "A4" and "A45" but one time the "A3". On the way back i used the combination "A4" and "A45" 2 times and the "A3" 3 times. The result should look like this:


texts:
I would like to simplify some polygons that represent rivers into linestrings. I assume that holes in polygons can be omitted. Anybody has some good idea how to simply do it?

It would be nice, if it was possible to be done directly in postgis or by another open source tool, but brief algorithm description will be sufficient.


texts:
I am able to connect to the mssql server database and choose the table that contains the spatial data. However, as soon as the layer is added, I get the 'crash dumped' error with a message saying
"minidump written to C:"

The database is hosted on azure
On a smaller database, I can see the nodes/location display briefly before the crash occurs.
I have tested the issue on a much smaller database that contained 1 row/entry and the crash still occurs but it's usually after I add another layer (e.g. open streetmap layer)


texts:
I am no expert in QGIS and I've already spend quite some time to find an answer to this, but I could not.

Could anyone tell me how to change the elevation value of DEM?

The one I have has been created for a small area, done with drone footage, and it has for some reason values that lie approx. 30 meters lower than the elevation a.s.l. really is (according to reference maps).


texts: Where to change GRASS GIS wxGUI language?

texts: I use GRASS 6.4.2 on Russian version of Windows 7 x64. As I start wxGUI it is shown in Russian. But I don't like its translation and want to work with English interface. I have no clues how to solve it.


texts: Qgis network analysis why find vertex return always -1

texts:
I'm trying to compute the shortest path between points using the QGIS networkanalysis library. But the graph.findVertex always return -1.

Basically what I'm trying to do is :
1/ Given a set of points, create a grid. The size of the grid is define by the poins bounding box
2/ transform all the cells' grid to point (centroid)
3/ For each cell's centroid compute the shortest path with the set of points

My script works fine with no error, but the graph.findVertex() function always return -1. It's like my points are not tied to the graph...

My code is available on pastebin :
http://pastebin.com/u6n4tSkz

Do you know what I'm doing wrong ?


texts: How do I delete empty fields from a shapefile with QGIS?

texts:
This question has been asked before but I couldn't get the solution to work. I have a shapefile with lots of empty columns and I would like to remove them all with a script or some tools.

This is the script that was posted as the solution to the other thread.

lyr = iface.activeLayer()
lstDelete = []

for idx in lyr.dataProvider().attributeIndexes():
    uv = lyr.dataProvider().uniqueValues( idx )
    if not uv:
        lstDelete.append( idx )

lyr.dataProvider().deleteAttributes( lstDelete )

# Update the layer structure
lyr.updateFields()

When I try it both on Windows 2.18.1 and Mac 2.18.0 I get no change to my active layer.


texts: How to convert EMF files to shapefile at QGIS? Actually I have EMF file that has been generated from Autocad-Civil3d.


texts:
I could not find anything on this topic already discussed.
We converted our file-based system in ESRI to an SDE instance on Sql Server.  We made changes in our code to open using the connection string and we change ArcMap to use the SDE as the default database.

We cannot open our GP services to republish them which is an indicator that ArcMap cannot find the schemas mentioned in the GP service ParameterInfo.  When I look at the database through ConfigurationManager the names have been changed from {schemaname} to {dbname}.DBO.{schemname}
and ESRI does not let me change the names back.

So do we have to alter every program to use the new schema name?  And are we locked into a given database name?  (If our client published to a different named server I can see us losing our ability to find the schemas once again)

Personally I think it is shortsighted to put the server name into the name of all of the fields in the ESRI database.  Surely there is a way around this???


texts: Blank spaces in large format plotter prints

texts:
I keep getting a blank swath when printing large format pdf files to my plotter.
I am trying to plot a single page that is 42" (maximum width of roll for this plotter) by up to 120" in length.  The map includes vector (lines, points, labels) on an image basemap (Bing Maps Road).
The pdf prints fine for 58" then there's about 4" of gap, then the plotter resumes.  I have experienced this on numerous files made from ArcGIS Pro and made from ArcGIS Desktop 10.2.1.  The pdf files look fine on the screen.  It seems to NOT be related to corrupt files or the writing of the pdf file, but rather problems in reading the pdf or sending it to the plotter.
I create the pdf in ArcGIS Pro 2.5.1 (or ArcMap 10.2.1).  I have tried deflate and adaptive compression.
I have tried 150 and 300 dpi.  I 'compress vector graphics', Max quality, embed fonts, no georeference, PDF layers only, RGB.
I avoid adobe reader.  I use Adobe Acrobat 11.0.23.
I have tried printing at 'Actual size" and at "100% of actual size"
I have tried normal quality, and best quality ('more passes' and NOT 'more passes').  I have tried printing as image (150dpi and 300dpi) and NOT as image.
I send to HP designjet T1700ps (driver 61.194.11903.100).  The driver is (@ Devices//Printer Properties/Advanced/) HP DesignJet T1700ps.  (Others in that list include PCL3GUI and HPGL2).
Please suggest a solution that will allow me to print plotter pages in excess of 58".


texts:
I have a gdb with survey data regarding urban trees. Each tree has photo attachments that correspond to the specific tree.

Below is a link to a ESRI script that will pull the photos and put them into a folder, but will not have them in any specific order, nor will they have a meaningful name. I would like to have a workflow that allows me to pull the attachments, rename them to something meaningful such as "456_344_1.jpg". This way I can easily navigate to a specific tree to confirm that it was identified correctly by the people who completed the survey.  The GDB has a unique id for each tree in its attribute table.

The goal is to populate a folder with photos of all of our trees. The format for file names should be something like: treeID_segmentID_uniquephoto.jpg

ESRI help article here: http://support.esri.com/technical-article/000011912


texts: How to display(expression) specific multiline labels (8 labels only) with same field in QGIS?

texts:
I wanted to display the labels (which only 8 labels with the longer names)from one field. I've tried in multiple ways but wasn't works (NOT duplicated question)

I've tried this:
if  ("NAMEASCII"='Nikolajewsk na Armur', 'Nikolajewsk-\nna-Armur',NAMEASCII)

[I've tried CASE statement, in my case works better only with IF statement]

The result of this expression (green: is the result, red: I need to do the same for this as well):

Layer properties > Label > Rule:

Label with:

BUT this is only one of them, there have 8 labels that I need to create the expression together (i don't know how to combine these 8 names only with one expression) just after the expression(above)
Can anyone tell me how? I may need a correct expression... because qgis giving me a syntax error, obviously, the expression that I created is wrong.


texts: Measuring distance between point to polyline using ArcMap

texts:
There are two lines

A : polyline
B : Axis

I created points along the line B with same interval using 'generated points along lines' tool of ArcMap.
I would like to measure the distance between these points and line A.
But If I use 'Near' tool, I could not get the perpendicular distance between line A and B.
How I can measure the ortho-distance between line A and points?


texts: How to LOLA DEM 16bit signed TIFF to QGIS Greyscale (BumpMap) image?

texts:
I have a moon 16bit signed TIFF downloaded to my computer from http://astrogeology.usgs.gov/search/details/Moon/LMMP/LOLA-derived/Lunar_LRO_LOLA_ClrShade_Global_128ppd_v04/cub

A pointer please in the right direction to convert the image to greyscale with QGIS but keep the greyscale in keeping with the altitudes.
That is to say running smoothly through the greyscale running white to black smoothly. Getting gradually darker as the altitude decreases?

In the downloaded TIFF...

WHITE is highest and then RED to YELLOW and on to GREENS, BLUES and VIOLETS.
Using Photoshop and choosing make Greyscale results in mixing the altitudes up.

Original picture source is huge.

I have downloaded 550 pixels per meter giving a file size of 1.43 GB.

Should be way more than i need to make a good bump map from the correctly converted greyscale.


texts:
I installed mapguide open source 2.6 and mapguide maestro 6.3 on windows 7 !! the tutorial that i followed assure that on this version of mapguide you don't need to install any server web. now when i try to go to : localhost:8008/mapguide it gives me an error of connexion, and the same error when i try with maestro. I note that the adresse : localhost:8008 work and i have apache folder and php folder in the directory where i installed mapguide..it's an  error on mapguide server i guess.

Can anyone help me?


texts: Extracting shapely polygons from .png file?

texts:
I have a .png file coming from a WMS service (WFS is not available). It only  contains "straight" lines and I would like to extract the geometries (shapely polygons) of the connected components of objects in the image.

I tried some basic solutions like using shapes from rasterio.features but the results were bad. There are 35 polygons in this image (if you pad the borders) and I either got thousands or a dozen. The polygons do not need to be perfectly aligned with the lines but I couldn't get even close to a working solution.

I consulted a solution to a similar problem (How to polygonize raster to shapely polygons) but it did not work well here.

I also tried finding contours and lines using various OpenCV methods but I still wasn't satisfied with the outcome. It seems like a pretty straightforward task but I can't find any working solution.


texts: Life cycle of an AxControls - license missing in some point of the execution

texts:
DEAR PEOPLE FROM THE FUTURE: Here's what we've figured out so far…

Does anyone know the lifecycle of an AxControl (AxTOCControl)?

I have a complex application, that uses several technologies (Spring.NET, NHibernate, DevExpress components) among them Esri's ArcObjects.

this application starts normally, it binds a product (engine or desktop) then it gets the licence (ArcEngine), and it starts, so it is getting the proper license.

Then, minutes after the aplication has started, it throws a dialog box saying that I do not have the license of ArcView (in fact I do not have it, but i do not need it, because my code only use functionalities related with ArcEngine, not even extensions).

If i click in the button of this dialog box, then the application is closed.

I have been trying to capture the thread that launches this dialog box, but this attempt has been useless.  Also I have tried using the tool Process Monitor and Process Explorer, but it has not been able to get information.

I need to understand the life cycle of the AxTOCControl to indetify when it is losing the license. Or if someone knows why a control can loss the binding to the existing license it will be very hepful!


texts: Changing symbology for layer using ArcPy?

texts:
I have a shapefile that is a 10-by-10 grid of squares.  The shapefile has 50 attribute columns (named p_1 through p_50) where each column has 100 entries (corresponding to the 100 grid squares) and each entry is one of 1, 2, 3, 4, or 5.

In order to visualize p_1, I navigate to Layer Properties, click on the Symbology tab, and click on 'Unique values' which is listed under Categories.  I then select 'p_1' from the Value Field, click the 'Add All Values' button, and click the 'Apply' button.  This results in a grid where each square has a color based on the number corresponding to the grid square in the p_1 column of the shapefile attribute table.  I then use the Export Map feature to save the colored grid as a png.

However, because there 50 columns corresponding to 50 different plans, I would like to automate this process.

How do I automate this process?

With help from @Fezter the following code accomplishes most of what I want:

import arcpy

mxd = arcpy.mapping.MapDocument("current")

# Find the plan layer
kdp_lyr = None
for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.name == 'kd_plans':
        kdp_lyr = lyr
        break

num_plans = 50
for i in range(0, num_plans):
    if kdp_lyr.symbologyType == "UNIQUE_VALUES":
        kdp_lyr.symbology.valueField = "all_the_plans$.plan_" + str(i)
        kdp_lyr.symbology.addAllValues()

        ## Remove the Null value
        sourceList = kdp_lyr.symbology.classValues
        sourceList.remove("<Null>")
        kdp_lyr.symbology.classValues = sourceList
        kdp_lyr.symbology.showOtherValues = False

    arcpy.RefreshActiveView()
    arcpy.mapping.ExportToPNG(mxd, 'plan_' + str(i) + '.png')
del mxd

The one thing it does not do is use the same color for each value (e.g., blue for the number 1 value across all plans i).

@Michael Miles-Stimson suggested using Apply Symbology from Layer however I don't believe this will work in my case because according to the documentation:

  The field in the Input Layer that will be displayed must have the same name as that of the corresponding Symbology Layer field. If this field is missing, the output data is drawn with default symbology.

and each field in my Input Layer has a different name.


texts: 78% cease calculation and "no message" from QGIS

texts:
I'm a total noob. I'm trying to generate an outline of a bunch of elevation lines which i can use to cut into some other elevation lines. I understand that I can buffer and dissolve to get the result I'm after, but when ever i try this the calculation runs to 78% and then stops and i get a very short popup message on my task bar that kindly says "no message."
I'm running an I7 12 core, 64 bit, 32GB ram, QGIS v 2.6.1.

Any suggestions would be welcome.


texts: Finding nearest point for each distinct category in performant way using big PostGIS table?

texts:
I have a table, etablissements with millions of rows, and a geom (st_point) column , spatially indexed (using gist).
Those features may have a different category, activite_principale.

I want to find from a given, fixed point the closest companies for 5 or 6 different categories called "activite_principale" (1 closest company per category).

Here's what I did right now:

    (WITH closest_candidates AS (
      SELECT
        ent.id,
        ent.name,
        ent.geom
      FROM
        geo_data.etablissements ent
      WHERE ent.activite_principale = '1071C'
      ORDER BY
        ent.geom <->
        'SRID=4326;POINT (5.4153978921979125 43.271437384501965)'::geometry
      LIMIT 10
    )
    SELECT id
    FROM closest_candidates
    ORDER BY
      ST_Distance(
        geom,
        'SRID=4326;POINT (5.4153978921979125 43.271437384501965)'::geometry
        )
    LIMIT 1)
    UNION ALL
    (WITH closest_candidates AS (
      SELECT
        ent.id,
        ent.name,
        ent.geom
      FROM
        geo_data.etablissements ent
      WHERE ent.activite_principale = '4711D'
      ORDER BY
        ent.geom <->
        'SRID=4326;POINT (5.4153978921979125 43.271437384501965)'::geometry
      LIMIT 10
    )
    SELECT id
    FROM closest_candidates
    ORDER BY
      ST_Distance(
        geom,
        'SRID=4326;POINT (5.4153978921979125 43.271437384501965)'::geometry
        )
    LIMIT 1)

--   UNION ALL
-- [...]  And so on...

I then clustered the etablissements table around the geom spatial index and ran VACUUM ANALYZE geo_data.etablissements;

Here's the result of EXPLAIN ANALYZE after clustering.

The planning is much shorter and the execution too but it's still slow (350-450ms).

I investigated compound index on geom & another text column, but that does not seem to be possible today ?

I use postgres 10 & postgis 2.4.

I don't know how to improve based on those explanations from the query planner. Can I do better performance ?


texts: Invalid Topology [Topoengine error] from ArcGIS for Desktop?

texts:
Using ArcGIS 10.1 for Desktop, I'm trying to run the Identity function on two large data sets but it fails every time. I don't think it is a lack of hardware as I can run other processes on the same data sets.

I've already updated software, and added the 64 bit geoprocessing add-in, checked and repaired geometry, and done the fix in KB29559 (set ARCTMPDIR env variable to a writable local folder).

Executing: Identity DudleyStamp OA2011 "D:\Case Studies\National2\National.gdb\DudleyStampOA11" ALL "1 Meters" NO_RELATIONSHIPS

Start Time: Wed May 08 21:11:43 2013
Reading Features...
Processing Tiles...
Assembling Tile Features...
ERROR 999999: Error executing function.
The table was not found.
The table was not found. [DudleyStampOA11]
The table was not found.
The table was not found. [DudleyStampOA11]
Invalid Topology [Topoengine error.]
Failed to execute (Identity).
Failed at Wed May 08 21:45:41 2013 (Elapsed Time: 33 minutes 58 seconds)

Based on comments so far I tried copying data to new file geodatabase, and running the tool from ArcCatalog still get exactly the same error.

I have ordered 16GB of the fastest RAM my computer can take (I current have 4GB). I also have a 6 Core AMD Phenom II X6 1090T, SSD, and AMD Radeon HD 6970 Graphics Card. So I'm not sure what more hardware I can add.

I've also tried increasing the size of the page file which was a suggestion I found on Google with no effect.


texts: ee.FeatureCollection.distinct() not working on certain properties

texts:
I did a batch download of the National Hydrography Dataset, and I have reason to believe that many lakes are found in many state datasets due to extending past state boundaries. I want to filter these out to avoid duplicates, but it seems to not work.
Below is a comparison of the sizes:
var nhd = ee.FeatureCollection("users/mincej20/NHD_Filtered_50");
nhd = nhd.sort("AreSqKm", False);

print(nhd.size()) //Returns 104762
print(nhd.size().distinct('NHDPlID')) //Returns 104762

These return the same sizes of featurecollections. However, if I print out the first 20 NHDPlIDs:
nhd = ee.FeatureCollection(nhd.toList(20, 1000));
print(nhd.aggregate_array("NHDPlID");

The following list is returned:
0: 22000600050688
1: 65000300111257
2: 65000300111257
3: 65000300111257
4: 65000200109644
5: 21000300292859
6: 21000100338800
7: 21000100338800
8: 21000500141259
9: 65000300137354
10: 65000300137354
11: 65000300137354
12: 65000200154805
13: 23000300045900
14: 23000300045900
15: 10000200330648
16: 15000300130262
17: 15000300130262
18: 15001700186644
19: 15001700186644
I'm not sure why this doesn't find distinct NHDPlIDs. I figured that I could find distinct '.geo' properties, and this seems to work, but I'd still like to figure out why finding distinct NHDPlIDs does not work. Here's some information about the features:
print(nhd.getInfo().columns);

AreSqKm: Float
Centroid: Geometry
Elevatn: Object
FCode: Integer
FDate: String
FType: Integer
GNIS_ID: String
GNIS_Nm: String
Max_Pixels: Long
Median_Occurrence: Float
NHDPlID: Float
Prmnn_I: String
ReachCd: String
Resoltn: Integer
Shap_Ar: Float
Shp_Lng: Float
VPUID: String
VsbltyF: Integer
system:index: String


texts: Unable to convert 3 band Spacenet dataset geotiff image to either jpg or png with gdal

texts:
I have 3 band images beloging to the Spacenet dataset which I'm trying to convert to jpeg. I get the following warning:

Warning 6: JPEG driver doesn't support data type UInt16. Only eight bit byte bands supported.

I then changed the ouput to png, but no matter what I do, gdal_translate outputs completely black images (checked with Gimp, it's all 0s). The geographic information seems to be extracted correctly.

gdalinfo provides the following info:

Driver: GTiff/GeoTIFF
Files: RGB-PanSharpen_AOI_2_Vegas_img1.tif
Size is 650, 650
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-115.307517599999997,36.128297699699999)
Pixel Size = (0.000002700000000,-0.000002700000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-115.3075176,  36.1282977) (115d18'27.06"W, 36d 7'41.87"N)
Lower Left  (-115.3075176,  36.1265427) (115d18'27.06"W, 36d 7'35.55"N)
Upper Right (-115.3057626,  36.1282977) (115d18'20.75"W, 36d 7'41.87"N)
Lower Right (-115.3057626,  36.1265427) (115d18'20.75"W, 36d 7'35.55"N)
Center      (-115.3066401,  36.1274202) (115d18'23.90"W, 36d 7'38.71"N)
Band 1 Block=650x2 Type=UInt16, ColorInterp=Gray
Band 2 Block=650x2 Type=UInt16, ColorInterp=Undefined
Band 3 Block=650x2 Type=UInt16, ColorInterp=Undefined

I've also tried with mogrify, but it doen's work either. In fact, it seems to get stuck.

I'm not sure whether I'm missing an intermediate step of some kind or there's something wrong with my file.

The image is here.


texts: How to concatenate double with string in field calculator with Python expression?

texts:
I'd like to concatenate my lat/lng with a few string fields in ArcGIS. What is the correct expression to use, avoiding VBA.  I don't want to convert the lat/lng fields to string either.

my 4 fields are category, type, x, y. with x and y being doubles.


texts: Getting Error: exported bands must be compatible from Google Earth Engine

texts:
I'm trying to export a tiff map and gave the following error:

  Error: Exported bands must have compatible data types; found
  inconsistent types: Float64 and Float32.

How can I solve this?

Code:

//Choose country using GEE Feature Collection

var region = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw').filterMetadata('Country', 'equals', 'Portugal');

//Add region outline to layer ‐ for selected countries

Map.addLayer(region,{}, 'Portugal');

// image collection pre 11/10/2017

var lt8_pre = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
  .filterBounds(region)
  .filterDate(ee.Date("2017-08-01"),ee.Date("2017-10-11"));

Map.addLayer(lt8_pre, {bands: ['B4', 'B3', 'B2'],min: 1000, max: 1500}, 'image L8 pre');

var lt8_ndvi_pre = lt8_pre
  .map(function(img){
  return img.addBands(img.normalizedDifference(['B5', 'B4'])).updateMask(img.select(['pixel_qa']).neq(5).neq(3).neq(2));
  });

// Create an NBR image using bands the nir and red bands (6 and 4)
var lt8_nbr_pre = lt8_pre
  .map(function(img){
  return img.addBands(img.normalizedDifference(['B7', 'B4'])).updateMask(img.select(['pixel_qa']).neq(5).neq(3).neq(2));
  });

print(lt8_ndvi_pre);
print(lt8_nbr_pre);

// image collection pos 11/10/2017

var lt8_pos = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
  .filterBounds(region)
  .filterDate(ee.Date("2017-10-18"),ee.Date("2018-02-10"));

Map.addLayer(lt8_pos, {bands: ['B4', 'B3', 'B2'], max: 1000}, 'image L8 pos');

var lt8_ndvi_pos = lt8_pos
  .map(function(img){
  return img.addBands(img.normalizedDifference(['B5', 'B4'])).updateMask(img.select(['pixel_qa']).neq(5).neq(3).neq(2));
  });

// Create an NBR image using bands the nir and red bands (6 and 4)
var lt8_nbr_pos = lt8_pos
  .map(function(img){
  return img.addBands(img.normalizedDifference(['B7', 'B4'])).updateMask(img.select(['pixel_qa']).neq(5).neq(3).neq(2));
  });

print(lt8_ndvi_pos);
print(lt8_nbr_pos);

//Clip to Specified Region
var NDVI_pre = lt8_ndvi_pre.mean().clip(region);
var NDVI_pos = lt8_ndvi_pos.mean().clip(region);
var NBR_pre = lt8_nbr_pre.mean().clip(region);
var NBR_pos = lt8_nbr_pos.mean().clip(region);

Map.centerObject(region, 10);
var ndvi_viz = {min:-0.8, max:1, palette:'000000,00FF00'};
Map.addLayer(NDVI_pre.select('nd'), ndvi_viz, "LT8 NDVI mean pre 11/10");
Map.addLayer(NDVI_pos.select('nd'), ndvi_viz, "LT8 NDVI mean pos 11/10");

// Display the NBRpost
var nbr_viz = {min: -1, max:1, palette: ['FFFFFF','CC9966','CC9900','996600', '33CC00', '009900','006600','000000']};
Map.addLayer(NBR_pre.select('nd'), nbr_viz, "LT8 NBR mean pre 11/10");
Map.addLayer(NBR_pos.select('nd'), nbr_viz, "LT8 NBR mean pos 11/10");

//Calcule dNBR and dNDVI
var dNBR = NBR_pre.subtract(NBR_pos);
var dNDVI = NDVI_pre.subtract(NDVI_pos);

//Map.addLayer(dNBR,nbr_viz, "dNBR");
Map.addLayer(dNDVI.select('nd'), ndvi_viz, "LT8 dNDVI");
Map.addLayer(dNBR.select('nd'), nbr_viz, "LT8 dNBR");

Export.image.toDrive({
      image: NBR_pos,
      description: 'NBR_pos',
      scale: 30,
      region: Pt,
      maxPixels: 278466375,
    });


texts: How can I store GPS and direction coordinates in EXIF data when geotagging photos on the Samsung S7

texts:
I am walking transects in the field for a research project and in addition to other GIS points I am loading photos as points into ArcMap using the Geotagged photos to points tool.

I have found that a Galaxy S7 smartphone is a great tool for taking the photos.
I notice however that the EXIF data from the photos I capture on my Samsung galaxy S7 has the GPS coordinates but that the direction field is empty.
I know that the phone measures this information and that Google Maps makes use of it.

Does anyone any experience of forcing the phone to store both GPS and Direction when geotagging ?


texts:
in short: My TIMESAT image-output of the seasonal parameter “Large Seasonal Integral” is wrong. TIMESAT continues to write a file with negative values (as 16 byte Integer) although Parameters should be positive and are positive in txt output files (TSM_printseasons). I do have the feeling that this is due to a scaling of the image or wrong dimensions in the header.

The “Small Seasonal Integral” image output is actually wrong as well, but if I define the image as UNSIGNED LONG INTEGER instead of a signed 16 Byte Integer as it supposed to be, the values of the output image are correct. This does not work with “Large Seasonal Integral”. Here , values are positive, but do not correlate to the txt output of TSM_printseasons.

This is the GUI view of my settings file within TIMESAT (looking quite good and neat ) :

This is the txt output of my time series which makes absolute sense. Please note the bigger Large Seasonal Integrals and the smaller Small Seasonal Integral (highlighted row is the 7th season Pixel 999.999 or 998,998 depending on which programm is used):

So if I open all Seasonal Parameters in ENVI, everything fits perfectly to the text output except the Integrals:

I smelled that something is wrong with the scaling, so I changed the header info to unsigned Long integer and then, the Small Integral was right as well. But not the Large Scale, although it changed to positive!

Now if I go back to TIMESAT and display (TSM_imageview) the binary image (meaning no header involved or whatsoever), I have to realize, that the binary image itself contains negative values! Although the text output of TIMESAT gives positive values!

So my colleagues and me are really puzzled. What did I do wrong? Why is it, that only one value is wrong? because the numbers are so large? How should the header look like? Might the problem sit there? But why is the binary image displayed in TIMESAT wrong itself?

Let me also add a screenshot of the header although I find it strange that the results WITHIN TIMESAT are different!


texts: Installing Python module within code and running from ArcGIS Pro script tool?

texts:
I am trying to install a third-party Python package from within a Python script so that other users can run the script and the package can be installed and loaded without having to manually download and install.

This is the code I'm using:

import subprocess

def install(package):
    subprocess.call([sys.executable, "-m", 'pip', 'install', package])

install('pyodbc')
import pyodbc

This seems to work fine from within my stand-alone Python script. However, when I try to run it as a script tool set up in ArcGIS Pro, it is not finding the installed package and is failing with this error message:

  ModuleNotFoundError: No module named 'pyodbc'

Do I need to change a path or environment setting?


texts: Dojo TypeError for ArcGIS Javascript API App

texts:
I recently updated the app to the 3.13 framework, and now am getting the following error when I load my application.

Uncaught TypeError: Cannot read property 'missingMessage' of undefined

Chrome says the error is being thrown on the Select Dijit.  I do not have "lang" in my HTML tag (like this error).  Anyone have any ideas?

var querySelect = new Select({
    id: "querySelect",
    name: "querySelect",
    options: [
        {label: "Signs", value: "signs", selected: true},
        {label: "Bridges", value: "bridges"},
        {label: "Culverts", value: "culverts"},
        {label: "Drains", value: "drains"},
        {label: "Roads", value: "roads"}
    ]
}, "query_select").startup();

UPDATE

Here is the simplifed code:

<script src="query.js"></script>    //Script with Select Dijit
<script>
    var map;

    require([
        "dojo/parser",
        "esri/config",
        "esri/IdentityManager",
        "dijit/layout/BorderContainer",
        "dijit/layout/ContentPane",
        "esri/map",
        "esri/toolbars/draw",
        "dojo/dom-construct",
        "esri/dijit/HomeButton",
        "esri/dijit/Scalebar",
        "esri/geometry/webMercatorUtils",
        "esri/dijit/BasemapToggle",
        "esri/dijit/Geocoder",
        "esri/layers/GraphicsLayer",
        "esri/layers/FeatureLayer",
        "esri/layers/ArcGISDynamicMapServiceLayer",
        "esri/Color",
        "esri/symbols/SimpleMarkerSymbol",
        "esri/symbols/SimpleLineSymbol",
        "esri/graphic",
        "esri/tasks/query",
        "esri/geometry/Extent",
        "esri/geometry/Point",
        "esri/tasks/GeometryService",
        "dojo/on",
        "dojo/dom",
        "dojo/dom-style",
        "dijit/form/CheckBox",
        "esri/tasks/IdentifyTask",
        "esri/tasks/IdentifyParameters",
        "esri/InfoTemplate",
        "esri/dijit/Legend",
        "esri/request",
        "dojox/grid/DataGrid",
        "dojo/store/Memory",
        "dojo/data/ItemFileWriteStore",
        "dojo/dom-construct"
    ], function (parser, esriConfig, IdentityManager, BorderContainer, ContentPane, Map, Draw, domConstruct, HomeButton, Scalebar, webMercatorUtils, BasemapToggle, Geocoder, GraphicsLayer, FeatureLayer, ArcGISDynamicMapServiceLayer, Color, SimpleMarkerSymbol, SimpleLineSymbol, Graphic, Query, Extent, Point, GeometryService, on, dom, domStyle, CheckBox, IdentifyTask, IdentifyParameters, InfoTemplate, Legend, esriRequest, DataGrid, Memory, ItemFileWriteStore, domConstruct) {

        parser.parse();


texts: QGIS 3x apply one data attribute table for more than 1 layer

texts:
I have got a few layers, in which the data attribute table looks exactly the same including:

Type (string)
Address (string)
Postcode (string)
X (East) (decimal)
Y (north (decimal)
Latitude (decimal)
Longitude (decimal)
ObjectID (string)
length (meters) (decimal)

I have to create them every time for every single layer. Is there any way to make them loaded instantly as a blank data attribute table for the layer I want?


texts: Importing CSV (downloaded from Google Drive) with polygonal geometry data, into QGIS

texts:
I formerly used Google Fusion Tables but due to the product downturn I had to bulk-download a backup copy of my data. Right now I have them all saved in CSV format.
I have attribute columns in these tables but also a column simply named as "geometry" all bunched up together, in a single column.
It's in the format of <Polygon><outerBoundaryIs><LinearRing><coordinates>[coord],[coord],[coord] [coord],[coord],[coord] [coord],[coord],[coord]</coordinates></LinearRing></outerBoundaryIs></Polygon> where [coord] is something like -129.83 for instance.
I'm trying to import this data into QGIS 3 and I've tried it both as a vector layer as well as delimited text file format. However, nothing actually shows up on the map.
What can I do to resolve this situation?


texts:
I would like to add some simple attribute validation in my current workflow. This is with regards to a point shapefile:

Required fields
Fields that become required if one optional field is populated
An 'ok' field that notifies the user that all required fields have
been populated

How would this be done?


texts:
I am trying to count raster values (from a .tif file) in polygons. In previous posts that people suggested using the 'zonal statistics' plugin in QGIS. However, I couldn't find it in my plugin manager (QGIS 3.6). I have ticked the box to show experimental plugins, but still couldn't find it. Any idea what is causing this?

If the plugin doesn't work, what would be a good alternative to use?


texts:
I would like to implement leaflet in shiny with R. I try to go with the simple examples I found online:

library(leaflet)
library(ggmap)

somePlace <- ggmap::geocode("Vienna")
# Information from URL : http://maps.googleapis.com/maps/api/geocode/json?       address=Vienna&sensor=false

somePlace
#   lon      lat
# 1 16.37382 48.20817

leaflet(somePlace) %>% addTiles() %>% addMarkers()

Map shows grey space with a blue marker. No trace of OSM spatial visualisation.

Could you point out where is the problem?

Similar problems (but without a suitable solution):

https://stackoverflow.com/questions/31070637/r-leaflet-package-produces-blank-grey-map-ubuntu-14-04
https://stackoverflow.com/questions/31647872/r-leaflet-rstudio-viewer-blank-grey-map-windows-7

Edit> sessionInfo() added on request:

R version 3.2.4 (2016-03-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows Server >= 2012 x64 (build 9200)

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252        LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
[5] LC_TIME=German_Germany.1252

attached base packages:
[1] grDevices utils     datasets  stats     graphics  methods   base

other attached packages:
 [1] rgdal_1.1-7          OpenStreetMap_0.3.2  osmar_1.1-7              geosphere_1.5-1      RCurl_1.95-4.8       bitops_1.0-6
 [7] XML_3.98-1.4         mapview_1.0.0        devtools_1.10.0      ggmap_2.6.1              leaflet_1.0.1        car_2.1-2
[13] zoo_1.7-12           dygraphs_0.8         plotly_3.4.1             dplyr_0.4.3          shinydashboard_0.5.1 shiny_0.13.2
[19] sp_1.2-2             R2HTML_2.3.1         Hmisc_3.17-2             ggplot2_2.1.0        Formula_1.2-1        survival_2.38-3
[25] lattice_0.20-33

loaded via a namespace (and not attached):
 [1] nlme_3.1-126        pbkrtest_0.4-6      xts_0.9-7           satellite_0.2.0     RColorBrewer_1.1-2  httr_1.1.0
 [7] tools_3.2.4         R6_2.1.2            rpart_4.1-10        DBI_0.3.1           lazyeval_0.1.10     mgcv_1.8-12
[13] colorspace_1.2-6    nnet_7.3-12         raster_2.5-2        gridExtra_2.2.1     quantreg_5.21       SparseM_1.7
[19] scales_0.4.0        hexbin_1.27.1       stringr_1.0.0       digest_0.6.9        foreign_0.8-66      minqa_1.2.4
[25] R.utils_2.2.0       base64enc_0.1-3     jpeg_0.1-8          htmltools_0.3.5     lme4_1.1-11         maps_3.1.0
[31] htmlwidgets_0.6     jsonlite_0.9.19     acepack_1.3-3.3     R.oo_1.20.0         magrittr_1.5        Matrix_1.2-4
[37] Rcpp_0.12.4         munsell_0.4.3       proto_0.3-10        viridis_0.3.4       R.methodsS3_1.7.1   stringi_1.0-1      [43] yaml_2.1.13           MASS_7.3-45         RJSONIO_1.3-0       plyr_1.8.3          grid_3.2.4              parallel_3.2.4
[49] rasterVis_0.37      splines_3.2.4       mapproj_1.2-4       rjson_0.2.15        gdalUtils_2.0.1.7   reshape2_1.4.1
[55] codetools_0.2-14    stats4_3.2.4        latticeExtra_0.6-28 png_0.1-7           nloptr_1.0.4        httpuv_1.3.3
[61] foreach_1.4.3       RgoogleMaps_1.2.0.7 MatrixModels_0.4-1  gtable_0.2.0        tidyr_0.4.1         assertthat_0.1
[67] mime_0.4            xtable_1.8-2        rJava_0.9-8         iterators_1.0.8     memoise_1.0.0       cluster_2.0.3  R version 3.2.4 (2016-03-10)


texts: Add spatial csv-layer in QGIS3 using Python3 code

texts:
I need to add a spatial csv-file with a specific path (folder and file name) in QGIS3x using a python script.
I have tried:
PyQGIS import csv- QGIS 3.4

uri = "file:///C:/users/heubacha/Test.csv?type=csv&delimiter=;&detectTypes=yes&xField=X_Wert_Pos&yField=Y_Wert_Pos&crs=epsg:25832&spatialIndex=no&subsetIndex=no&watchFile=no"
mem_layer = QgsVectorLayer(uri, 'Test', 'delimited text')
QgsProject.instance().addMapLayer(mem_layer)

no luck with that...warning: invalid data provider delimited text.

Trouble loading non-spatial data via CSV using Python for standalone script
was not getting me anywhere either, because I need to add the path as a string somewhere in the coding, and that's not described there.


texts: Calculate length of Multilinestrings/Linestrings/Polygons in geoJSON file

texts:
I have a couple of routes on my map where I would like to display the length of each route. I think I've found what I was looking for here: https://github.com/tyrasd/geojson-length, but I'm not sure how to apply the calculations.

I've tried to refer to the files I've downloaded with:

<script src="js/measureLengthGeoJSON.js"></script>

but then I get this error message:

Uncaught ReferenceError: module is not defined
    at measureLengthGeoJSON.js:1

Line 1 refers to: module.exports.

I then tried to copy the entire formula/function like so:

var lengthGeo = function(geometry) {
    if (geometry.type === 'LineString')
        return calculateLength(geometry.coordinates);
    else if (geometry.type === 'MultiLineString')
        return geometry.coordinates.reduce(function(memo, coordinates) {
            return memo + calculateLength(coordinates);
        }, 0);
    else
        return null;
}

function calculateLength(lineString) {
    if (lineString.length<2)
        return 0;
    var result = 0;
    for (var i=1; i<lineString.length; i++)
        result += distance(lineString[i-1][0],lineString[i-1][1],
                           lineString[i  ][0],lineString[i  ][1]);
    return result;
}

/**
 * Calculate the approximate distance between two coordinates (lat/lon)
 *
 * © Chris Veness, MIT-licensed,
 * http://www.movable-type.co.uk/scripts/latlong.html#equirectangular
 */
function distance(λ1,φ1,λ2,φ2) {
    var R = 6371000;
    Δλ = (λ2 - λ1) * Math.PI / 180;
    φ1 = φ1 * Math.PI / 180;
    φ2 = φ2 * Math.PI / 180;
    var x = Δλ * Math.cos((φ1+φ2)/2);
    var y = (φ2-φ1);
    var d = Math.sqrt(x*x + y*y);
    return R * d;
};

and apply the function like so:

var testlength = lengthGeo(HvenKort);

console.log(testlength);

Nothing is outputted to the console. I can't figure out what I'm doing wrong.

Any suggestion on how to fix it or any other library that is capable of calculation distances of Multilinestrings/Linestrings/Polygons?


texts: Problem in refreshing FeatureLayer in Openscales

texts:
I have a refresh feature layer problem, because I need to delete all the markers and add the new markers so I tried this piece of code, but it doesn't work :(

            //refresh
            markers.clear()
            markers.reset();
            markers.removeFeatures(markers.features);
            markers.redraw(true);
            //refresh

            trace(markers.features,"features before");
            trace(Snaps ,"Snaps" );
            for each(var Car in Snaps)
                {
                    var mat = Car[0];
                    var Loc = Car[2];
                    var snap_id = Car[1];
                    var direction = Car[3];
                    //Alert.show(Loc[0]);
                    var marker:CustomMarker = CustomMarker.createUrlBasedMarker("images/markers/iconFlag1up_.png", new Location(Number(Loc[1]),Number(Loc[0])));
                    marker.addEventListener(MouseEvent.CLICK,handlesClickMarker);
                    marker.name=snap_id;

                    markers.addFeature(marker);//add marker to markers
                }

                trace(markers.features,"features after add");


texts: TIN interpolation errors and weird lines [Qgis 2.18]

texts:
I want to create an elevation model in EPSG 4326 from 2 shapefiles [img0] : elevation lines [Altimetria linee Antrodoco] and elevation points [Altimetria punti Antrodoco], both obtained by the District Maps database of Lazio.

However I encounter several errors and/or weird rasters as a result, depending on my Interpolation settings:

a) if I choose both Points and Lines as Input [img 1 in the link above] and set them as Structure lines and Points respectively, I get minidump error [img 2], no matters the cellsize or columns values.

b) if I choose both Points and Lines as Input [img 3], and set both of them as Points, I do not get errors, but I get alwais messy/weird lines, that do not seem to follow the elevation lines at all [img4]. This happens also with clipped shapefiles.

c) all in all, whenever I choose points or lines, as soon as I change the cellsize to get 5 (project unit is meters), the rank and columns change to 0. [img 5] Only with values inferior to 0.00x I get positive ranks and columns value [img 6], independently from the size of the layer.

Since I'm a beginner, I believe there are probably many errors and issues braided together. What is the first/main problem in my procedure?

I can't add more than two links, so I put all the images in the first gdrive folder.


texts: Error loading netCDF data using sf package in R

texts:
I'm trying to use netCDF data, and would like to use sf rather than ncdf4 if it's already built to load netCDF data into sf objects.

However, I'm finding that reading the file using sf::st_read isn't working.

You csn get the data for the below example here: ftp://ftp.cdc.noaa.gov/Datasets/cpc_global_temp/tmax.2018.nc

df <- st_read("tmax.2018.nc",)
# Reading layer `tmax.2018' from data source `tmax.2018.nc' using  driver `netCDF'
# Warning message:
# no simple feature geometries present: returning a data.frame or tbl_df

It can't find any of the geometries or the data. I'm guessing this is because I am not specifying the layer correctly.

Any help as to how to find the right layer or how to get sf to play nicely with netCDF data?


texts: Difference between Positioning and Localization

texts:
I would like to know what are the differences between Positioning and Localization. In most papers they are used interchangeably.

Are they the same?
Does one come after another?

I always thought in this order:

Localization: GPS coordinates found.
Positioning: Found you in a map.
Navigation: Now you can navigate thru a map.

Is this correct?


texts: How to insert a list of parameters from a parametric view in OpenLayers?

texts:
After successfully (with the precious help of members from stackexchange) creating a Geoserver parametric view that works as it should, now i'm facing another problem:
How to pass a list of parameters from a given view in a new OpenLayers.Layer.WMS call.
From the geoserver user manual it is stated that The viewparams argument is a list of key:value pairs, separated by semicolons:

viewparams=p1:v1;p2:v2;...

But since i'm calling the viewparams in an openlayers new wms layer, as soon as i put the semicolon on the script:

 var x = new OpenLayers.Layer.WMS(
        "y",
            "http://10.194.32.103:8080/geoserver/wms",
            {'layers': 'vista', 'format':'image/png', 'transparent':'true', viewparams:'parameter:data';'parameter1:data1'},
            {'opacity': 1, 'isBaseLayer': false, 'visibility': false}
            );

i get an error:

Any ideas on how to overcome this situation?

Thanks!


texts: Creating point attribute distribution map in QGIS

texts:
I'm working with multiple species distribution data and am looking at creating what amount to density distribution maps of those species. I have a small number of site locations (18) and only 10 species that are being tracked. Each site has the presence of from one to 10 species in values ranging from 1 - 59.

I want to create a layer that shows the density distribution of all of the species at the same time. Heatmaps in QGIS only do point density (is my understanding at least, maybe I'm wrong) and this can be weighted by an attribute, but isn't very accurate. I've done an IDW interpolation and then created contour lines, but the IDW algorithm looses the data from the low value sites and then creating contours excludes the low value sites (i.e. sites with just 1 of the species) as the values are not evenly distributed. See below.

Maybe I'm overcomplicating this, but is there a way to show the distribution of point attribute values that is more accurate and doesn't involve an interpolation that essentially erases the lower values?


texts:
I've been trying to struggle through Terrset but just can't get it to do what I want. I initially had problems downloading the files from geoexplorer but seem to have figured that out. The problem now is that I deleted my original project file with the intention of starting fresh. However, I'm now getting errors whenever I open Terrset and also when I try to start a new project. When opening Terrset I get the following errors:

Data path not found. Please select another project or change the data path.

There is an inconsistency in the user's preference file "C:\Program Files (x86)\Terrset\TerrSetUsers.ini", line "projectdir=C:\Users\jdleh\Documents\Romote Sensing\lab2pt2". That folder does not exist

When I try to select a new project in the project tab I get the following error:

Cannot create file "C:\Users\jdleh\Documents\Remote Sensing\lab2pt2\lab2pt4.env". The system cannot find the path specified

Any thoughts on these errors?


texts: Why does iface.activeLayer() not work in a PyQGis standalone app?

texts:
It is my second day on PyQGIS. Not sure why layer doesn't return anything. Can anyone help me on this?

Code:

#!/usr/bin/python
from qgis.core import *
from qgis.utils import *
from PyQt4.QtCore import *

def main():
  qgis_prefix="/usr"
  QgsApplication.setPrefixPath(qgis_prefix, True)
  QgsApplication.initQgis()
  app = QgsApplication([], False)

  uri = "elevp.csv?delimiter=%s&xField=%s&yField=%s&elevField=%s" % (";","x","y","elev")
  vlayer = QgsVectorLayer(uri, "test", "delimitedtext")
  QgsMapLayerRegistry.instance().addMapLayer(vlayer)
  #print(QgsMapLayerRegistry.instance().mapLayers())

  #iface.mapCanvas().setSelectionColor( QColor("red") )
  layer = iface.activeLayer()
  print(layer)

# Run main wrapper codes
main()


texts: Convert Leaflet WMS layer to Leaflet GeoJSON layer

texts:
Using the Leaflet API I am adding a WMS layer through the L.tilelayer.WMS feature. This feature includes only a certain amount of options such as 'layers','crs','style' etc.

Through the Leaflet documentation I figured that the L.GeoJSON feature includes the 'onEachFeature' option. I would like to use this 'onEachFeature''  option within the WMS layer.

What would be a solution to apply an 'onEachFeature' option within a Leaflet WMS TileLayer?


texts: Removing layers from OpenLayers not working after layers are drawn

texts:
I have a simple map application written with Angular and Node.js
There is search service that works fine and shows the results with geometry.
After the search results I draw each result on map
So far everything is fine.
But after each result, I have to remove the old layers drawn on the map which I can not do it.
Here is my code:
onSearch(params: any) {
   if (params.target.value !== "") {
       this.httpClient.get('http://127.0.0.1:3001/search?q=' + params.target.value).subscribe((response: any) => {
       this.searchResults = response
       this.mapService._map.removeLayer(this.vectorPolygon) // removing layers
       this.drawLocation(this.searchResults) // draw new ones
       return response
    })
   } else {
       this.searchResults = []
    }
}
drawLocation(results: any) {
results.map(res => {
  const format = new ol_format()
  const feature = format.readFeature(res._source.st_astext, {
    dataProjection: 'EPSG:4326',
    featureProjection: 'EPSG:3857'
  })
  this.vectorSource = new ol_source_vector({
    features: [feature]
  })
  this.vectorPolygon = new ol_layer_vector({
    source: this.vectorSource,
    displayInLayerSwitcher: false,
    style: this.styleMe(feature)
  })
  this.mapService._map.addLayer(this.vectorPolygon)
})
}


texts: Renaming field in attribute table via PyQGIS?

texts:
How can I rename a field in PyQGIS?

I tried:

for field in layer.pendingFields():
    if field.name() == 'old':
        layer.startEditing()
        oldname = field.name()
        field.setName('new')
        newname = field.name()
        print(u'Old name: {}, New name: {}'.format(oldname,newname))
        layer.commitChanges()

This code print 'Old name: old, New name: new', but in attribute table is field  name still 'old'.


texts:
I have been trying for a couple of days to get hillshading to work as I need it to but have run out of things to try and have tried the solutions I have found through searching the web.

I've created a TIFF and imported into QGIS and applied any number of combinations of hillshading settings but the result always ends up looking like the image on the left. The right handside shows the resolution and level of detail of my TIFF (tried 72dps upto 600dpi) and still the same. Is there a setting I am missing that will enable me to get a smoother gradient with my hill shading?

(I should add the 2 images are not showing the same area of the map, they are from the same map, just not the same area. It was a screen quick shot)


texts: QGIS: Joining BCGS .csv color HEX file to BCGS polygon .shp

texts:
I know that topic has been brought up, but I still don't get it.
Please, see the 2 pictures below: what do I missed?
It seems like I can join the .csv but the colors still don't show correctly. The only colors working are from the QGIS categorized random ramp. I want the standard BCGS colors from the .csv to be applied to the .shp.
Thanks for your help,


texts:
I am unable to get a ResultImageLayer.

I can get a GPFeatureRecordSetLayer using ResultData and add the feature geometries individually.

I am using ArcServer 10.1 so it is not a pre-release problem.

I can see the map service that accompanies the geoprocessing service and the legend / symbology in this map service.

I am doing this with the JavaScriptAPI.

Using Firefox I see a “500 Internal Server Error” but there is no other messages.
Could be as simple as my call:

gp2.getResultImageLayer(jobInfo.jobId, "Test01", imageParameters, addGPServiceResult2, errorProcessing)
function addGPServiceResult2(gpImageLayer) {
   map.addLayer(gpImageLayer);
}

I’m not sure what the ParameterName “Test01” should refer to although I’ve tried ever possibility that seemed reasonable.

Or could be an ArcServer setting?

Any help will be greatly appreciated.

Thanks,
Charles


texts: QGIS, creating layer from PostGIS : add a style

texts:
I am creating some layers in QGIS using PostgreSQL connection.

I have a python plugin doing the load of the layer in Qgis,
this part of the code looks like this :

        uri = QgsDataSourceURI()
        db = SettingsProperties.get_db()
        pg = SettingsProperties.get_postgre_settings()
        uri.setConnection(db['host'], str(db['port']), project_name, pg['user'], pg['password'])
        uri.setDataSource('work', table_name, str(column_geom))

        layer = QgsVectorLayer(uri.uri(), tr(layer_name), "postgres")

No problem with this, but now, I would like to setup style when creating this layer, this relying on the value of some column ... etc

Actually, I am able to setup style through a .qml prepared,
but I wanted to know if there is a possibility to create the style either in PostgreSQL or in the Python script in a "simple" way and then add it to the database ?

I know you can save the style with QGIS GUI. But I would like to script it.


texts: Finding the best locations for bee hives within a specified area (QGIS)

texts:
I would like to find fixed positions for bee hives within an orchard to ensure complete pollination of the trees. Is there a way to do this in QGIS?

Example: The orchard area is 2 hectares (I have a polygon dataset representing the orchard). 10 bee hives are required to pollinate the area. How would I go about calculating the best positions for the the bee hives to ensure complete bee coverage of the orchards in QGIS?


texts: How to remove small areas without reducing borders

texts:
I'm having a problem where I am overlaying two maps to create a new kind of map with the data conjoined. In general this is good, but I need to remove the smaller areas from it. For the interior polygons, you can see that this is fine. However, for the ones on the exterior, they are removed and area is actually lost.

I don't want this to happen - I need to persist the original borders. I would like extended shape to have the data values of one of the polygons it was previously connected to.

Original map:

Cleaned map:

Additional side note... There's an issue with v.buffer in Grass 6.4.3 where it will get in some sort of an infinite loop, use 100% of a CPU core, and never complete. I would like to avoid doing any operations that involve buffering the areas of the original map. However, if that's the only way please let me know. I will try to work around this problem.

Thanks!


texts:
I need to set cell values in one raster to NoData based on another raster. Both rasters have the same extent and some cells overlap. For those cells that overlap I want to set them to NoData in Raster1 while keeping the rest of the cell values.

Please see graphic: Raster1 (green); Raster2 (black). I need to replace the cell values in Raster1 that overlap Raster2 with NoData while keeping the rest of the cell values in Raster1.

I'm using the SetNull function in Raster Calculator but I'm getting a blank raster so I guess the expression is not quite right.

SetNull("%Maas_Vol_47_com%","%Maas_Vol_47_inc%")
Raster1 = Maas_Vol_47_com
Raster2 = Maas_Vol_47_inc


texts: Georeferencing 1990 U.S. county names / FIPS for CartoDB?

texts:
I have a data table with columns of county and state names as well as FIPS codes for 1990 U.S. counties.

I cannot easily merge shape data into the data (e.g. the Census ones), but CartoDB does not recognize the "administrative regions" being these counties at all when I ask to georeference this data.

Is there a way around this?


texts:
How can I set the zoom to 23 or more? When I try to set the zoom higher the map image disappears. I am able to zoom in past 23 once the map is activated in a browser, but I can't set the zoom past 23 to begin with. Thoughts?
var map = L.map('mapdiv', {zoomControl: true}).setView([40.712,-74.006], 22);

L.tileLayer('https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}', {
  maxNativeZoom:18,
  minZoom: 4,
  maxZoom: 35,
}).addTo(map);

L.control.scale().addTo(map);


texts: How show filter image out of feature in other color?

texts:
I am using OpenLayers JS to filter the features display. I am trying to show the filtered imaged in other color keeping rest of feature also.

Following is the code I am using for filtering:

var a="73.84";
var b="73.89"
var filter= new OpenLayers.Filter.Logical({
     type: OpenLayers.Filter.Logical.AND,
   filters: [
    new OpenLayers.Filter.Comparison({
       type: OpenLayers.Filter.Comparison.BETWEEN,
         property: "x",
         lowerBoundary:a,
 upperBoundary:b
            }),
            new OpenLayers.Filter.Comparison({
                type: OpenLayers.Filter.Comparison.BETWEEN,
                property: "y",
                lowerBoundary:"18.51",
                upperBoundary:"18.56"
            })]
    });

var parser = new OpenLayers.Format.Filter.v1_1_0();
var filterAsXml = parser.write(filter);
var xml = new OpenLayers.Format.XML();
var filterAsString = xml.write(filterAsXml);

allwards.params["FILTER"] = filterAsString;
allwards.redraw();

Using this code my feature is filtered but original feature are removed and only filtered feature are show. I want to show original features along with filtered feature in different color.

Please help me in this case.


texts: VectorTileServer works using web links but not using ArcMap--Add ArcGIS server

texts:
In ArcMap I am trying to add the following ArcGIS server
https://staging.gcm.telstra.com.au/gisserver/rest/services/Hosted/Mobile_3G_ExtAntenna/VectorTileServer
It connects fine but I don't see anything in the hosted folder

If I use a link from the page above I can see the service in a web map
https://staging.gcm.telstra.com.au/gisserver/rest/services/Hosted/Mobile_3G_ExtAntenna/VectorTileServer?f=jsapi

Is there any way to see this in ArcMap (10.8)


texts: Use of the Peclet number in FEFlow as a mesh-quality metric

texts:
The user manual for FEFlow (v7.0) states

  "In transport models, the elemental distribution Péclet Number is
  available which can be used to check whether the spatial
  discretization is sufficiently fine for a specific transport problem."

However, it is not clear whether a higher or lower Peclet number indicates an appropriate mesh for the transport problem.

Can anyone shed some light on this?


texts: How can I change mouse cursor while it hovers over the map in Oracle Mapviewer

texts:
I'm using Oracle Mapviewer. I want to switch mouse pointer from default to crosshair when user enters "Interactive mode". I'm developing in GWT and tried do do something like this:

public static void setCursorCrosshair() {
    removeCursorStyles();
    getMapWidget().addStyleName("cursorCrosshair");
}

But it changes mouse pointer only in few parts of map, it looks like it omits map tiles and change cursor to crosshair on FOI's that I've added just before.
How can I make it to display crosshair over the whole map?


texts: Subtracting out whole regions using occurrences (points)

texts: I have a bunch of occurrence points of animals scattered across a global map, and an eco-region map underneath it. I want to identify the eco-regions (i.e. Tundra) that does NOT contain any occurrence points. How can I go about doing that?


texts: What is a good resource for making choropleth maps using excel data?

texts:
I would like to make a choropleth map of Malaysian regions, colour-coding them by population density.

Is there a place where I can find templates of such a map?


texts: arcpy.da.SearchCursor ignores where_clause parameter?

texts:
I'm trying to use arcpy.da.SearchCursor to go through records in a feature class stored in a SQL SDE database.  I do so through a feature layer (created in arcpy) called "unit_lyr".  When I run the da.SearchCursor, the where clause is ignored and arcpy returns all records in the layer.  When I run it using the old arcpy.SearchCursor, everything runs as expected.

Any idea why?

# Run it the old way, everything goes fine.
unit_list = []
rows = arcpy.SearchCursor('unit_lyr', 'UnitDate IS NOT NULL')
for row in rows:
    if row.UnitName is not None:
        if len(row.UnitName) > 1:
            unit_list.append(row.UnitName)

>>> len(unit_list)
192

# Use arcpy.da and all values are returned.
unit_list = []
with arcpy.da.SearchCursor('unit_lyr', 'UnitName', 'UnitDate IS NOT NULL') as rows:
    for row in rows:
        if row[0] is not None:
            if len(row[0]) > 1:
                unit_list.append(row[0])

>>> len(unit_list)
1085


texts: Average of multiple rasters in PostGIS raster?

texts:
I am trying to calculate the average of multiple rasters stored as raster datatype in PostGIS enabled database. The rasters are stored in three different tables (spi3, sri3, and smdi) and have date as common column. I am trying to use ST_Union function of PostGIS, where it says we can specify setof raster.

This is what I have done so far.

select st_union(SETOF[spi3_raster, sri3_raster, smdi_raster], 'MEAN')
from
(
select spi3.raster as spi3_raster, sri3.raster as sri3_raster, smdi.raster as smdi_raster
from spi3 as spi3
join sri3 as sri3
  on spi3.date = sri3.date
  and sri3.date = '2013-01-01'
join smdi as smdi
  on sri3.date = smdi.date
where spi3.date = '2013-01-01'
) as foo;

However, I am not getting the result as expected.


texts: I'm really confused by this. ArcMap will seeminly randomly decide to not allow snapping when editing. Sometimes its only for certain features, sometimes for only parts of the same feature (e.g. one end of a line but not the other, or not along a part of the line). Now its just flat out not letting me snap to anything (yes, snapping is turned on with all options on). Any hints on this one?


texts: add a shapefile layer to a specific postgis table

texts:
I have a question about Postgis:
There is a way(extension) to import a shapefile and add it directly to a specific table (not database), if not there is a request to joint vertically two tables like : i have table A and table B

                                 Table "A"
         gid             |             Name             |        geom
-------------------------+------------------------------+-------------------
          1              | Arganier                     | xxxxxxxxxxxxxx

                                    Table "B"
         gid             |             Name             |        geom
-------------------------+------------------------------+-------------------
          1              | Abousier                     | xxxxxxxxxxxxxx

And then i want to joint this two tables and get the following table:

                                 Table "A"
         gid             |             Name             |        geom
-------------------------+------------------------------+-------------------
          1              | Arganier                     | xxxxxxxxxxxxxx
-------------------------+------------------------------+-------------------
          2              | Abousier                     | xxxxxxxxxxxxxx


texts:
I want to import several dxf-features. This works fine in QGIS 2.18 apart from text (annotations).

The respective text information is stored, at least in my case, within the point components of dxf features and can be displayed in QGIS by utilizing labels.

I would prefer if the labels would be displayed directly (like in CAD software and even ArcMap) since angles of text do not seem to be stored.

Is there a workaround for showing dxf text components automatically and/or for visualizing text labels at correct angles? (dxf2shp did not help me either)


texts: How do I create an outline layer for multiple blocks of points?

texts:
I have a series of points that mark the edges of excavation units, I would like to create a layer of an outline connecting the points. The end goal is a layer displaying a grid of the excavation units. I've tried the points-to-line feature, but the only success is creating horizontal lines using Northing and Easting as line and sorting features. At one point, I was able to use the selection and construction tools to complete the outline and have it selected, but not only could I not save that data as a layer or new shapefile, I have four separate blocks of excavation units marked and starting on another block automatically cleared the selected grid.

I'm including a photo showing the situation, the red Xs mark the corners of units, they are each a single point.

I would like a layer that is a grid, demonstrated in part by the black lines drawn on the picture.


texts: IdentifyTask for "hollow" polygon (using ArcGIS API for JavaScript)

texts:
I have a web map application set up that's currently testing on a variety of Dynamic Map Service layer types, including points, lines, and polygons. This application includes a tool that makes use of ItentifyTask to allow a user to click on the map and get a popup with attribute information for any visible layers at that point. This tool works perfectly for points, lines, and SOLID polygons, but it does NOT work for polygons that are displayed with a hollow symbology (ex: county polygons displayed to only show the outline of the county). In this case, the execution of the IdentifyTask returns NO results. If I click on the outline of the feature it works fine, but not if I click on the hollow interior.

Is there a good workaround for this? For our needs, we really do need the ability to identify such layers, and the hollow symbology is our standard for such layers.


texts: Grass gis - 3D data set interpolate to 3D raster

texts:
I have a dataset of 70,000 points and I would like to interpolate the properties as a 3D raster.  I have created 'slices' of properties each meter of elevation in the model domain; approx: 200 slices.

The format of the dataset is as follows:
BH, EastingNad83, NorthingNad83, Elevation, Property(numerber),Property(text)

(coordinates are in UTM)
The data points are located along boreholes (approx 5000).

I am able to import the data set in Grass Gis. However, I am unable to create the 3D raster using V.Vol.Rst.  Is the data set too large for this module?


texts: QGIS doesn't render raster alpha channel properly (jagged edges)

texts: I have two rasters in my project clipped by custom mask and translated to jpg. One is perfectly fine, but the other one has jagged edges which appear right after translating & saving to jpg. Photoshop sees both rasters correctly, but QGIS renders only one properly, while the other one has those saw-like edges. I tried other source rasters, other masks, other QGIS versions and nothing works, besides saving to tiff. Does anyone know how to fix this?


texts: Fill in no data polygons using data from surrounding polygons

texts:
I'm using ArcGIS 10.4 and I have an advanced license.

I have a polygon shapefile that has 100000 polygons. The attribute table has two columns Id and soil type.

In my data, there are 1000 polygons with NO attributes regarding soil type.

How can I use the common soil type in the nearby polygons to fill in the the attributes in polygons with no soil type?


texts: Prevent gdalinfo from creating auxiliary file

texts:
When calling gdalinfo with certain options, it creates an auxiliary file which has a side-effect on other options.
Let’s say I want to compute the min/max value on elevation data from a GeoTIFF file. I am using the -mm option:
gdalinfo /path/to/file.tif -mm

It returns the usual gdalinfo output with image size, projection, etc. plus the following due to the -mm option:
  Computed Min/Max=592.000,2737.412

So far so good.
I can also compute approximate stats:
gdalinfo /path/to/file.tif -approx_stats

Which adds this output:
  Minimum=611.662, Maximum=2737.412, Mean=1461.410, StdDev=414.722
  Metadata:
    STATISTICS_MAXIMUM=2737.412109375
    STATISTICS_MEAN=1461.4103933815
    STATISTICS_MINIMUM=611.66223144531
    STATISTICS_STDDEV=414.72166382554

Values are slightly different, probably due to the approximation, but it’s not an issue for me.
However, this call also writes an auxiliary file next to the original file, e.g. /path/to/file.tif.aux.xml
This has several downsides. First, it creates a file on a folder where I should not be writing into. Secondly, it now pollutes basically every call to gdalinfo on that same file.
For instance if I call gdalinfo with -mm exactly like the first call, the statistics are always output, mixing the approximate and "real" min/max:
  Min=611.662 Max=2737.412   Computed Min/Max=592.000,2737.412
  Minimum=611.662, Maximum=2737.412, Mean=1461.410, StdDev=414.722
  Metadata:
    STATISTICS_MAXIMUM=2737.412109375
    STATISTICS_MEAN=1461.4103933815
    STATISTICS_MINIMUM=611.66223144531
    STATISTICS_STDDEV=414.72166382554

I am currently deleting the auxiliary file after every call to gdalinfo that writes one, but I would prefer a solution that does not write that file to begin with.
PS. adding the -nomd option (no metadata) suppresses metadata from the console output, but the auxiliary file is still written to disk.


texts: How to view values from stacked raster maps of a particular pixel at the same time (ArcGIS)

texts:
Re: ArcGIS v10.1

I have a stack of raster maps.  I would like to click on a pixel, and view ALL of the values for that particular pixel at the same time.  The only tool that I could find is the pixel inspector, which only displays the value for one map at a time.  I want to compare all values from my layered maps for this particular location at the same time.

Has anyone done this with ESRI software?


texts: GeoServer 2.10.0 - How to use CreateStoredQuery

texts:
I need to create the GetSpatialDataSet stored query in my WFS service to follow the INSPIRE rules.

That WFS service is in a GeoServer 2.10.0. I have checked the official docs about how to create a stored query: http://docs.geoserver.org/stable/en/user/services/wfs/reference.html

The problem is I don't know what I am supposed to do with that piece of code.

<wfs:CreateStoredQuery service='WFS' version='2.0.0'
 xmlns:wfs='http://www.opengis.net/wfs/2.0'
 xmlns:fes='http://www.opengis.org/fes/2.0'
 xmlns:gml='http://www.opengis.net/gml/3.2'
 xmlns:myns='http://www.someserver.com/myns'
 xmlns:topp='http://www.openplans.org/topp'>
<wfs:StoredQueryDefinition id='myStoredQuery'>
<wfs:Parameter name='AreaOfInterest' type='gml:Polygon'/>
<wfs:QueryExpressionText
 returnFeatureTypes='topp:states'
 language='urn:ogc:def:queryLanguage:OGC-WFS::WFS_QueryExpression'
 isPrivate='false'>
  <wfs:Query typeNames='topp:states'>
    <fes:Filter>
      <fes:Within>
        <fes:ValueReference>the_geom</fes:ValueReference>
         ${AreaOfInterest}
      </fes:Within>
    </fes:Filter>
  </wfs:Query>
</wfs:QueryExpressionText>

Where should I paste it? It says it's a POST request, but that doesn't tell me anything.

EDIT: I have tried the following script from another answer:

<html>
<head>
<title>WFS Test</title>
<script language="JavaScript">
  function doit()
  {
    var url = 'http://myserver:8080/geoserver/test/wfs';
    var method = 'POST';
    var postData = '<wfs:CreateStoredQuery\n'
    + 'xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"\n'
    + 'xsi:schemaLocation=\n'
    + '"http://www.opengis.net/wfs/2.0\n'
    + 'http://schemas.opengis.net/wfs/2.0/wfs.xsd"\n'
    + 'xmlns:cp_inspire="urn:x-
inspire:specification:gmlas:CadastralParcels:4.0"\n'
    + 'xmlns:fes="http://www.opengis.org/fes/2.0"\n'
    + 'xmlns:wfs="http://www.opengis.net/wfs/2.0"\n'
    + 'xmlns:gml="http://www.opengis.net/gml/3.2"\n'
    + 'xmlns:xsd="http://www.w3.org/2001/XMLSchema"\n'
    + 'service="WFS"\n'
    + 'version="2.0.0">\n'
    + '<wfs:StoredQueryDefinition\n'
    + 'id=\'' +
'http://inspire.ec.europa.eu/operation/download/GetSpatialDataSet' + '\'>\n'
    + '<wfs:Parameter name=\'' + 'CRS' + '\' type=\'' + 'xsd:string' + '\'/>\n'
    + '<wfs:Parameter name=\'' + 'DataSetIdCode' + '\' type=\'' + 'xsd:string' +'\'/>\n'
    + '<wfs:Parameter name=\'' + 'DataSetIdNamespace' + '\' type=\'' + 'xsd:string' +'\'/>\n'
    + '<wfs:Parameter name=\'' + 'Language' + '\' type=\'' + 'xsd:string' + '\'/>\n'
    + '<wfs:Parameter name="count" type="xsd:integer"></wfs:Parameter>\n'
    + '<wfs:QueryExpressionText\n'
    + 'returnFeatureTypes=\'' + 'cp_inspire:CP.CadastralParcel' + '\'\n'
    + 'language=\'' + 'urn:ogc:def:queryLanguage:OGC-WFS::WFS_QueryExpression' + '\'\n'
    + 'isPrivate=\'' + 'false' + '\'>\n'
    + '<wfs:Query typeNames=\'' + 'cp_inspire:CP.CadastralParcel' + '\' srsName="${CRS}">\n'
    + '</wfs:Query>\n'
    + '</wfs:QueryExpressionText>\n'
    + '</wfs:StoredQueryDefinition>\n'
    + '</wfs:CreateStoredQuery>\n';

var req = new XMLHttpRequest();
req.open("POST", url, true);
req.setRequestHeader('User-Agent', 'XMLHTTP/1.0');
req.setRequestHeader('Content-type', 'test/xml');
req.onreadystatechange = function () {
  if (req.readyState != 4) return;
  if (req.status != 200 && req.status != 304) {
    alert('HTTP error ' + req.status);
    return;
  }
  alert(req.responseText);
}
if (req.readyState == 4) return;
req.send(postData);

}

But I'm getting this error:

Does anyone see where is the problem?


texts: Deleting multiple vertices at once from representations using ArcGIS Desktop?

texts:
I have a polygon feature class that I've converted to representations so that I can delete some vertices.

I'm doing this to simplify some of the shapes manually, but I'd like to speed this up by selecting a group of vertices with the Representations lasso tool and then deleting them at once rather then using the time-consuming point-and-click Delete Vertex tool.

Is this possible?


texts: Create a rough "world overview" map with mkgmap

texts:
I'm trying to generate a "World overview map" for my Garmin device. The idea is to have this map on the device all the time, for a very rough idea of the position and then feed in higher resolution maps where necessary, e.g. when going on trip. I would like to use OSM for this purpose with mkgmap.

If I do this straightforward, I would create a ridiculously large img-file, using all the OSM data. I would thus like to just have lower scale information included in the map. However I haven't grasped yet, how to achieve this. Can anybody help me out?


texts: Preparing hillshade map in mapinfo discover

texts: Is there a tool to perform hillshade in MapInfo Discover like there is in ArcGIS Desktop?


texts: Looking for topographic map of Bogota, Colombia, with elevation contour lines

texts:
I'm looking for a topographic map of Bogota, Colombia; specifically, we want to get the elevation contour lines, and export that layer to CAD.

If anyone knows of any free places to find this data, we would love to hear it! I've been looking through some of the sources here http://freegisdata.rtwilson.com/ - but I'm not very familiar with that huge list.

Does anyone know of a specific link in there that would be best to investigate?


texts: Using two Python Addin Combo Boxes to choose Folder then Layer

texts:
I am trying to create a Python add-in for ArcMap for users to quickly select layer files from within several different folders. This will require 2 combo boxes. User selects the folder from the first combo box, then the layerfile from the second combo box. The second combo box is populated based on the folder selected in the first.

I have found a couple of questions on here similar to this: Using two Python Addin Combo Boxes to choose feature class first then field from chosen feature class? and Dynamically adding items into one python addin Combobox from another?

I have successfully copied and pasted PolyGeo's code for selecting feature class then field, but when I try and change the code to select folder then layer file, I can't get it to work. The first combo box displays the folders correctly, but the second combo box is blank. Can anyone help?

Code below. I am a beginner with ArcPy and add-ins.

import arcpy
import pythonaddins

class ComboBoxFolders(object):

    """Implementation for AddIns_addin.combobox (ComboBox)"""
    def __init__(self):
        arcpy.env.workspace = r"C:\\EditData\\Folders\\"
        self.items = arcpy.ListWorkspaces("*", "Folder")
        self.editable = True
        self.enabled = True
        self.dropdownWidth = 'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
        self.width = 'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
    def onSelChange(self, selection):
        ComboBoxLayers._hook.items = [x.name for x in arcpy.ListFiles(selection)]

class ComboBoxLayers(object):

    """Implementation for AddIns_addin.combobox (ComboBox)"""
    def __init__(self):
        self.items = []
        self.editable = True
        self.enabled = True
        self.dropdownWidth = 'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
        self.width = 'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
        ComboBoxLayers._hook=self
    def onSelChange(self, selection):
        mxd = arcpy.mapping.MapDocument("CURRENT")
        addLayer = arcpy.mapping.Layer(r"C:\\EditData\\Folders\\" + selection)
        df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
        arcpy.mapping.AddLayer(df, addLayer, "TOP")


texts:
I was looking for a way to add a new column to several shp (simple polygons) at the same time. I don't know much about batch command.

This new attribute is also a little bit tricky : I need to generate new datas which will be a classification based on the standard deviation from an existing column.


texts: How to set individual feature angle in OpenLayers

texts:
I have an application where each element has an angle. I can use "Angle" function in "style" but not individual. Is there any way to apply angle individually?

function DesenhaPonto(Id, Estilo, Longitude, Latitude, Source, Angulo) {
            try {

                var coordinates = [Longitude, Latitude];
                var feat = new ol.Feature(new ol.geom.Point(coordinates));
                feat.setStyle(Estilo);
               // feat.angle(Angulo); <---- We need set individual Angle here.
                feat.name = Id;
                Source.addFeature(feat);
            } catch (e) {
                TrataErro(e);
            }
        }


texts: How to understand/interpret slope unit values in QGIS 2.14?

texts:
When I create slope raster I get this result:

I can't understand the those values represent?

What are the units of those values?

I saw this question QGIS slope percent calculation , but didn't find answer.


texts: Maps with numbers and postcodes using QGIS?

texts:
I have a list of UK postcodes and a number of participants in a survey. I would like to map the first three letters of the UK postcode and I would like a shape....a circular shape for example which will show a figure and its size is relative to the figure, ie. 50% is half the size of 100% in each postcode area.
This will show how many people participated from different areas.

I am totally new to QGIS and I was hoping that I could get a link to either a guide on how to do this, or link to practice data together with a step-by-step guide.


texts: I have a tif file that needs some pixel values adjusted in ArcMap. I have a pixel with a value of 44 that I need to change to 58.  In the Raster Calculator I'm using Con("tif" == 44,58,"tif") but that changes the lowest value that I need to keep the same.  The high value can change though if need be.  Still pretty new to spatial analysis.


texts: What size-reduction to expect on average when gzipping MBtiles?

texts:
As part of checking some prereqs for an idea (IOS app with MBtiles on local storage enabeling offline maps) I'd like to know if it's useful to gzip MBtiles when transporting them to an IOS device.

In other words, is there a useful reduction in size when gzipping MBtiles? (or is the MBtile-format, already packed in some way, thus limiting the use of gzip or other packers). If so, how much size reduction can I expect? (percentage ballpark)

NOTE: crosspost from SO (as it's within both domains) https://stackoverflow.com/questions/13864691/what-size-reduction-can-i-expect-on-average-when-gzipping-mbtiles


texts: Why are certain GCPs in my georeferenced scanned map resulting in warping? (QGIS)

texts:
I am having issue with georeferencing a scanned map. Warping has occurred in the West in Morocco (near Marrakesh) as well as a sort of "black holes" in Libya and Northern Chad. As far as I can tell, the GCP coordinates are correct, and I have checked to see if the reference/scanned map is accurate, and as far as I can tell, it is. Furthermore according to documentation for QGIS (which is what I am using) the dX and dY fall within acceptable variation/range (not being greater than 10), which presumably rules out the GCPs being misplaced.

Any idea as to what may be the issue? The transformation setting for target SRS maybe? Other portions of the map are referenced fine.


texts: How to Access Secured ArcGIS Online Services from a Custom JS Web app?

texts:
I have some Layers Hosted on my company's ArcGIS Online Account. I want to make a Custom WebMapping Application hosted on our public server, which General Public can use without Logging in.

I'm trying to find some document or Example which shows how this can be done, but have drawn a Blank.

My Research so far:

1)Accessing secure ArcGIS Online map service from custom web mapping application? does not answer my question because I don't want my users to login to AGOL

2) https://developers.arcgis.com/documentation/core-concepts/security-and-authentication/accessing-arcgis-online-services/ Doesn't show how to create an App on an Organisation's ArcGIS Account.

What's the correct way to Access Secured ArcGIS Online Services from a Custom JS Web app?


texts: Dynamic Map Service with Zoom Slider in arcgis by API for javascript

texts: I have ArcGISDynamicMapServiceLayer and I want to show Zoom Slider by API 3.2 for javascript


texts: Creating automated tests using Google Earth Engine Python API

texts:
I'm creating a lib using the Google Earth Engine Python API.
I need to test this script in travis CI. When I'm launching the unittest into my env it work but when I push to Github the test fail because the virtual env doesn't know my Google credential.
Is there a way to add my Earth Engine credential into my test environment?


texts: Masking Water Bodies to build surface temperature raster slice of them in ENVI?

texts:
I am looking at a region in the Alaskan tundra, there are plenty of small water bodies in my scene.

I need to mask these bodies, so that I can perform a surface temperature/emissivity normalization scene on only these water bodies. The ultimate product is to look at the original scene, and overlay the temperature masks unto the bodies of water.

How can I do this?

Edit: I am relatively new to ENVI, but I have tried doing a NDWI with band math.. It created a scene in which it identified possible water bodies in white and the rest in black. I could not figure out how to create a shapefile out of the white bodies it identified..

I have also considered classifying the water bodies but I am not sure how to go about this and ultimately create a mask.

I am working with Landsat 8 Data.


texts: Making layertree with multiple nested nodes?

texts:
On Adding subnodes to layermanager? I asked how to get two levels on a layermanager. The response provided there was what I was needing at that moment. Now, I´ll like to create a little more complex structure in the layer tree. The expected structure is something like:

|overlays|
|        |layer1
|        |layer2
|Group1
|      |layer3
|      |layer4
|      |Subroup1
|      |        |layer5
|      |        |layer6
|      |
|      |Subroup2
|              |layer7
|              |layer8
|background
|          |layer9
|          |layer10

I found Adding multiple nested nodes in GeoExt? doing something similar but to be honest I don’t understand the code, and don’t know how to integrate it into my own code to use it with the suite-sdk.


texts: Weighted Distance with Focal Statistics / Raster Calculator

texts:
I have a raster image of a road network with pixel values corresponding to vehicule frequency.

What I would like is to calculate a weighted distance raster with the frequency as weight within a maximum distance. I don't seem to grasp the concept of how to reference neighbouring cell values.
How do I refine and apply a formula like the following or better in the raster calculator AND / OR how can I make Focal Statistics work with a weighting formula?
newraster(x,y) =
     for all i in range(-maxdist, maxdist):
        for all j(-maxdist to maxdist):
            oldraster(x-i,y-j)*(1-distance/maxdistance)


texts: QGIS export attribute table to CSV and keep value relation

texts:
If I export my joined table to CSV I get the numeric id value instead of the value relation name from the joined table.

Can I export my table with the value relation value?


texts: CRS conversion with epsg_operations.properties btw. EPSG:32041 and EPSG:3857

texts:
I try to correct the visualisation of EPSG:32041 on a Bing map. For that purpose, when i import the 32041 data i use a correction in geotools with epsg_operations.properties.

I debugged, which transformation is created when no entry for 32041,3857 is present in epsg_operations.properties and took that as a base which i than specified in epsg_operations.properties.

However, without changing any of the parameters in the transformation the data imported without an entry in epsg_operations.properties differs from data imported with entry in epsg_operations.properties.

I debugged and found out that the transformation in epsg_operations.properties is not used directly, but its further transformed. That happens in PropertyCoordinateOperationAuthorityFactory.createFromCoordinateReferenceSystemCodes(String, String, boolean) Small part of it here, mt is the transformation from properties file, mt2 is the transformation used for CRS reprojection:

    // Create MathTransform from WKT
    MathTransform mt = null;
    try {
        mt = factories.getMathTransformFactory().createFromWKT(WKT);
    } catch (FactoryException e) {
        // Probably malformed WKT.
        LOGGER.warning("Error creating transformation: " + WKT);
        return null;
    }

    // Create the CRS definitions
    String s = this.authority.getIdentifiers().iterator().next().getCode();
    CoordinateReferenceSystem source = CRS.decode(s + ":" + sourceCRS);
    CoordinateReferenceSystem target = CRS.decode(s + ":" + targetCRS);

    // Need to create a derived MathTransform that will handle axis order and units
    // as defined in CRS. Had to cast to DefaultMathTransformFactory because
    // createBaseToDerived is not defined in MathTransformFactory interface (GeoAPI).
    DefaultMathTransformFactory mtf =
            (DefaultMathTransformFactory) factories.getMathTransformFactory();
    MathTransform mt2 = mtf.createBaseToDerived(source, mt, target.getCoordinateSystem());

DefaultMathTransformFactory.createBaseToDerived() sais, it expects a "raw" transform without unit conversion or axis switch.
Maby this is my problem? If yes, how do i get this raw transformation in WKT?
Thanks a lot!


texts:
I would like there to be only 3 levels of zoom available in the WMS that is displayed in my OpenLayers.Map (from level 6 to level 8).

demolayer2 = new OpenLayers.Layer.WMS(
    "abc...","http://localhost:8080/geoserver/gwc/service/wms",
    {
      layers: 'def...',
      transparent:"true",
      format: 'image/png',
      numZoomLevels: 3,
      minZoomLevel: 6},
   {isBaseLayer:false},
   {tileSize: new OpenLayers.Size(256,256)}
);

map.addLayer(demolayer2);

map.addControl(new OpenLayers.Control.LayerSwitcher(
   {'div':OpenLayers.Util.getElement('layerswitcher')})
);

map.zoomToExtent(new OpenLayers.Bounds(
    -4.601615515076983,39.8769407866263,-3.0527184873294764,41.16710732525929)
);


texts: Merge two adjacent polygons in a 3D space


texts:
In Leaflet I would like to assign different icons to my Point of Interests based on attribute "tipo" in geojson file:

{ "type": "Feature", "properties": { "Name": "Fontana dei Colli", "tipo": "fontana", "img_id": "colli", "description": "Via Campasso", "image_url": "img/fontane/FontanaDeiColli.jpg", "timestamp": null, "begin": null, "end": null, "altitudeMode": null, "tessellate": -1, "extrude": 0, "visibility": -1, "drawOrder": null, "icon": null }, "geometry": { "type": "Point", "coordinates": [ 10.48693469784657, 45.39481875814014, 0.0 ] } },

This is my js code:

function PoIstile(feature, layer) {
switch(feature.properties["tipo"].toString()) {
        case "fontana":
            return {
            icon : L.icon({
                iconUrl: 'icons/fontana-small.png',
                iconSize:     [22, 30], // size of the icon
                iconAnchor:   [22, 22], // point of the icon which will correspond to marker's location
                popupAnchor:  [-3, -26] // point from which the popup should open relative to the iconAnchor
            }),
            }
            break;
        case "chiesa":
            return {
            icon : L.icon({
                iconUrl: 'icons/chiesa-small.png',
                iconSize:     [22, 30], // size of the icon
                iconAnchor:   [22, 22], // point of the icon which will correspond to marker's location
                popupAnchor:  [-3, -26] // point from which the popup should open relative to the iconAnchor
            }),
            }
            break;
            }
layer.bindPopup("
"+feature.properties.Name+"

");
};
var PoI = L.geoJson(PoI, {
onEachFeature: PoIstile
});

It doesn't work but it doesn't return any errors in console, but I can't find the problem


texts: Mapnik pgraster in PostGIS using ST_MapAlgebra on Tiled rasters

texts:
Is it possible and feasible to run a simple ST_MapAlgebra on the fly for a multi-band tiled (100x100) raster with mapnik extension pgraster?

For example I am trying with this xml:

<Layer name="raster-layer" status="on" srs="+proj=utm +zone=56 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs">
<StyleName>raster-style</StyleName>
<Datasource>
    <Parameter name="host">localhost</Parameter>
    <Parameter name="dbname">data</Parameter>
    <Parameter name="user">gisuser</Parameter>
    <Parameter name="password">pass</Parameter>
    <Parameter name="table">(select ST_MapAlgebra(rast, 1, rast, 2, '([rast1.val]+0.00000001)/([rast2.val]+0.00000001)') as rast from maps.my_raster as r where rast is not null) as foo</Parameter>
    <Parameter name="type">pgraster</Parameter>
    <Parameter name="raster_table">maps.my_raster</Parameter>
    <Parameter name="raster_field">rast</Parameter>
</Datasource>
</Layer>
</Map>

When I request a tile, the request takes several minutes or times out so I am guessing my request is looking at too much data. I was hoping I could get the request to run the map algebra on just the 100x100 pixels row before passing it on to pgraster.

I have looked at OptimizeRenderingWithPostGIS which helps but is directed more that the vector data types and I think the tile rasters adds another layer of complexity.

If I run the below query in psql it takes about 1 second for 1 100x100 tile.

select rid,ST_MapAlgebra(rast, 1, rast, 2, '([rast1.val]+0.00000001)/([rast2.val]+0.00000001)') as rast from maps.my_raster as r limit 1


texts: ModelBuilder: Error 999999 "Failed to copy raster dataset"

texts:
I want to run the "Project Raster" tool in the ModelBuilder of ArcMap but I get the "ERROR 999999: Error executing function. Failed to copy raster dataset."
and I have no reason why.

Input is a simple Tiff raster and the projection is looped like described here:
Looping Projection in ArcGIS ModelBuilder?


texts: Building CityGMLs with OSM, OSM2World and FME

texts:
Here is my use-case, I am looking to write a singular FME workflow that takes XMLs from Overpass (full of building ways, parts and relations), runs them through OSM2World to get OBJs and then finally converts them to CityGML with FME.

I have it working in principal (https://github.com/stirringhalo/osm2citygml) and with tweaks to the Overpass query, I have it working for either all buildings in a bounding box, or all building parts of a relation, etc.

My question is this: What is the best approach in FME to either split an OBJ such that each building is a feature (so that each bounding box isn't a single CityGML feature) or failing that, iterate every building using Overpass or similar so that I can feed it into OSM2World and FME? Or would it be simpler to not merge geometries but then how would I eventually add attributes to buildings if each feature is disjoint?

I'd like to avoid running each building through FME directly as I've noticed the startup time for FME is slow, so it would accumulate pretty dramatically for a large number of buildings

The FME workspace is here https://raw.githubusercontent.com/stirringhalo/osm2citygml/master/osm2citygml.fmw